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Abstract of “System Identification and Model Reduction Using Modulating Function Tech- 
niques” by Yan Shen, Ph.D., Brown University, May 1993 

Weighted least squares (WLS) and adaptive weighted least squares (AWLS) algorithms are 
initiated for continuous-time system identification using Fourier type modulating function 
techniques. Two stochastic signal models are examined using the mean square properties 
of the stochastic calculus: an equation error signal model with white noise residuals, and a 
more realistic white measurement noise signal model. The covariance matrices in each model 
are shown to be banded and sparse, and a joint likelihood cost function is developed which 
links the real and imaginary parts of the modulated quantities. The superior performance 
of above algorithms is demonstrated by comparing them with the LS/MFT and popular 
predicting error method (PEM) through 200 Monte Carlo simulations. A model reduction 
problem is formulated with the AWLS/MFT algorithm, and comparisons are made via six 
examples with a variety of model reduction techniques, including the well-known balanced 
realization method. Here the AWLS/MFT algorithm manifests higher accuracy in almost 
all cases, and exhibits its unique flexibility and versatility. Armed with this model reduction, 
the AWLS/MFT algorithm is extended into MIMO transfer function system identification 
problems. The impact due to the discrepancy in bandwidths and gains among subsystems 
is explored through five examples. Finally, as a comprehensive application, the stability 
derivatives of the longitudinal and lateral dynamics of an F-18 aircraft are identified using 
physical flight data provided by NASA. A pole-constrained SIMO and MIMO AWLS/MFT 
algorithm is devised and analyzed. Monte Carlo simulations illustrate its high-noise reject- 
ing properties. Utilizing the flight data, comparisons among different MFT algorithms are 
tabulated and the AWLS is found to be strongly favored in almost all facets. 
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Chapter 1 


Introduction 

1.1 Overview 

Mathematical models of dynamical systems are often required in engineering, physics, 
medicine, economics, ecology and in most areas of scientific enquiry. In control and 
system engineering, model-building or system identification from measurements of 
input-output data on a dynamical system has been one of the most active fields 
drawing enormous attention from researchers around the world. Among the many 
well established parameter estimation schemes, algorithms like the prediction error 
method (PEM) [1] [22] [12] have enjoyed a sustained boom in the past decade for 
discrete-time models. Although many researchers have been utilizing different kinds 
of transformations in an effort to link both continuous and discrete time model iden- 
tification into a single framework, grim problems persist, like nonuniqueness of the 
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transformations, making these methods unwieldy and potentially unreliable. Hence, 
a direct attack on the problem is clearly preferred if a continuous time differential 
model is desired. The research presented in this thesis will be focused on the develop- 
ment of methodology in parameter identification for linear continuous-time differential 
equation models. 

Generally, the identification of linear differential systems can be formulated in a 
deterministic vein using the classical steady state frequency domain approach for 
estimating the system transfer functions, or using a variety of time domain methods 
like Bellman-Kalaba’s quasilinearization [3] [16], state variable filters, model reference 
techniques and adaptive observers [56]. In a stochastic vein, the known methods 
would include generalized least squares, instrumental variables, maximum likelihood 
and extended Kalman filtering techniques [56] [1]. 

Stemming from Shinbrot’s method [51] [28] of moment functionals and using a set 
of carefully chosen modulating functions to facilitate converting a differential equa- 
tion on a finite time interval into a set of algebraic regression equations in param- 
eters, the modulating function technique (MFT) has been exploited as a tool for 
identifying continuous-time models. The modulating process itself can be viewed 
as discretizing a continuous-time differential system into a corresponding frequency 
domain characterization by means of a “resolving frequency” u Q . Pearson and Lee 
[37] [18] [38] [35] utilized a set of real commensurable sinusoids {cosmu;ot,sinmu;of}, 
m = 0, 1,2, • • ■ , M, where loq = 2-k/T is the resolving frequency or the “step size” in 
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the frequency domain, to build up the modulating function set through solving Van- 
dermonde type linear equations 1 with which differentiating the data and estimating 
unknown initial conditions for time limited data can be totally avoided. By contrast, 
the other forms of modulating function methods, e.g., Hermite polynomials used by 
Takaya [55], Poisson process by Fairman and Shen [7] and Saha and Rao [43] [44] 
[45] [46], require either a long time interval of data or constrained initial conditions. 
Computationally, the modulating process by Fourier type modulating functions can 
be efficiently implemented by well documented Fast Fourier Transformation (FFT) 
algorithms while the other algorithms have to face a heavier numerical burden asso- 
ciated with the process of converting a differential equation on a finite time interval 
into a set of algebraic equations. 

Besides the earlier work mentioned above, Co and Ydstie [6] have applied the trigono- 
metric based Fourier type modulating function technique (FTMFT) to model re- 
duction and some MIMO system identification problems 2 in chemical engineering. 
Pearson and Pan [36] [26] further expanded the FTMFT into the nonparametric 
identification framework, under which three least squares nonparametric algorithms 

for estimating system transfer functions are formed. Shen and Pearson [49] applied 

x For numerical reasons, singular value decomposition (SVD) is strongly recommended [49] for 
the higher order modulating function sets. Pearson has suggested a much more efficient complex 
form of modulating function set (the set used in this thesis) which totally avoids tackling these 
Vandermonde linear equations. 

2 The identified MIMO models presented in that paper have higher orders than the actual MIMO 
systems, i.e. , some unobservable modes were included in the models thus obtained. 
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the trigonometric-based MFT to analog Butterworth filter banks used in wind tunnel 
experiments in fluid dynamics, formulated Kalman filter type recursive schemes for 
both parametric and nonparametric LS/FTMFT algorithms, implemented the nu- 
merically sound Bierman’s U-D matrix factorization algorithm for updating the high 
dimension Kalman gain matrix, developed the Parzen- window- based order determina- 
tion algorithm for parametric LS/FTMFT, tested a parallel-adaptive memory-saving 
paradigm for the nonparametric LS/FTMFT method used in the on-line configura- 
tions, and thoroughly demonstrated the affect of higher order assumptions for both 
model and the modulating function set on the quality of the final estimation. Using 
the nonparametric approach, Pan [26] developed frequency matching model reduction 
algorithms which performed better than the parametric LS/FTMFT used by Co and 
Ydstie [6]. A FTMFT-based high resolution frequency estimation method [27] for 
signal processing has been proposed. Meanwhile the FTMFT has also been applied 
to the time-varying systems and nonlinear system identification problems [40] [29] 
[38] [31] [30] [32] [34] [33]. For brevity, FTMFT will be abbreviated as MFT in the 
rest of this thesis. 


1.2 Organization of This Thesis 

Following a brief introduction to MFT early in Chapter 2, one fundamental mod- 
ulating property will be established showing that a modulated time domain white 
Gaussian stationary process will be a stationary Gaussian stochastic sequence in the 


4 



discrete frequency domain with its covariance matrix being banded by the order of 
the modulating function set. Then an idealistic equation error signal model is in- 
troduced leading to the first weighted least squares (WLS/MFT) algorithm which is 
based on a maximum likelihood estimate. For the much more realistic stochastic sig- 
nal model with additive measurement noise, the explicit form of the regression error 
covariance matrix, which is a function of the unknown parameters, will likewise be 
shown to be banded but no longer stationary in the discrete frequency domain. Using 
this covariance matrix as a weighting, a numerical relaxation scheme dubbed as the 
adaptive weighted least squares (AWLS/MFT) algorithm will be devised, which is 
an approximated maximum likelihood estimate 3 . The third part of Chapter 2 deals 
with important implementation issues incurred by the use of a complex modulating 
function set. In order to combine both the real and imaginary parts into a unified 
cost function, one Lemma regarding vital relationships among the different covari- 
ance matrices is established as a prerequisite to constructing a joint likelihood cost 
function linking the information from the modulated real and imaginary quantities. 

In order to assuage the affliction caused by inverting the covariance matrix, a recur- 
3 The maximum likelihood estimate here is different from the ML estimate in [39] where multi- 
nonoverlaping data blocks are required, but in the AWLS/MFT, only one single data block is needed 
and this is perhaps a more economical and realistic framework especially in time limited transient I/O 
data. Meanwhile the framework of the AWLS/MFT can also guarantee a finer resolving frequency 
for the same short length of data. As a matter of fact, the AWLS/MFT algorithm was partly inspired 
and initiated so as to ease the urgency of efficiently utilizing rather limited available flight data in 
an aircraft identification problem (see Chapter 5). 
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sive banded-sparse matrix inversion algorithm will be derived and its computational 
efficiency and stability will be elucidated. Finally, a 200 Monte Carlo simulation 
comparison study is presented for the AWLS/MFT, WLS/MFT and LS/MFT, and 
a comparison is made with the popular prediction error method (PEM), where the 
improvement and superiority of WLS/MFT and AWLS/MFT algorithms will be illu- 
minated. This chapter is the most important part of the thesis and will serve as the 
theoretical cornerstone for the rest of the thesis. 

Armed with AWLS/MFT, Chapter 3 deals with model reduction. The performance 
of AWLS/MFT will be evaluated for six examples published in the literature pro- 
viding comparisons with other known methods like nonparametric frequency fitting 
[26], FF-Pade approximation [14], time-domain optimization [47] and the well known 
Balanced-Realization (B-R) technique [24]. As one of the two best model reduction 
schemes in our comparison studies, AWLS/MFT is found to be able to perform at 
least as well as the B-R scheme, and in addition, possesses a kind of flexibility and 
versatility that the B-R algorithm lacks. 

As the second application of the AWLS/MFT algorithm, the most general setting 
of MIMO system identification will be considered in Chapter 4. Inasmuch as the 
original transfer function matrix might not be modulatable directly in this case, one 
procedure is suggested for converting from an unmodulatable form to a higher order 
modulatable differential system. Based on the measurement noise signal model, the 
decomposability from MIMO into a set of MISO models is discussed from the view- 
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point of a joint likelihood cost function. AWLS/MFT is applied to each member of 
the set of MISO systems to get an identified model for the higher order system 4 , and 
then the model reduction scheme of Chapter 3 is employed to obtain the original un- 
modulatable transfer function matrix. Five MISO systems with different bandwidth 
and magnitude combinations are used as examples to illustrate the impacts of these 
combinations on the accuracy of estimations and the overall feasibility and applica- 
bility of this approach. Results from 200 Monte Carlo runs under moderate additive 
noise settings have been quite encouraging. 

In Chapter 5, flight data of an F-18 aircraft provided by NASA will be utilized, as 
a comprehensive application example of the WLS/MFT and AWLS/MFT algorithm, 
to identify the longitudinal and lateral dynamical systems of the aircraft. Due to the 
physical constraints posed by the aircraft itself, AWLS/MFT is first extended into a 
pole-constrained form, contrary to the decoupled form used in Chapter 4, and then 
employed to tackle the physical flight data. Simulation studies on these algorithms 
also manifest good noise- rejection features. Other extended AWLS/MFT algorithms 
based on the coupled state space models are devised as well, though they do not 
produce as good results as the I/O-based constrained AWLS/MFT algorithms. 


4 Co and Ydstie [6] used LS/MFT to accomplish this step and took it as the final result which 
actually leads to higher order models containing unobservable modes. 
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Chapter 2 


Weighted Least Squares in MFT 


2.1 Brief Introduction to the Basics of MFT 

Consider the following nth order SISO differential equation system 

“n 6 n _,u (,) (<) + e(t), a 0 = 1 (2.1) 

t =0 1=0 

where {a^} and {6,},z = 1,2, are the time-invariant parameters needed to be 

identified from the input-output data pair {it(f),y(f);t € [0, 7 1 ]} and superscript “(z)” 
means zth derivative, i.e., y^°\t) = y(t ) and y^(t) = d'y(i)/dt u , e(t) represents the 
effect of modeling errors. Assuming smoothness, the property which an n th order 
modulating function 4>(t ) has to satisfy relative to a fixed time interval [0, T] is: 

^W(O) = <f>W(T) = 0; i = 0, 1,2, ...,n — 1. (2.2) 
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Multiplying both sides of (2.1) with <f>(t) and then integrating by-parts over [0, T], 
while noting (2.2), leads to the following essential relation of the MFT: 

E(-l)‘ a n-i [ T y(i)<f> {i) (t)dt = £(-l yb n -i [ T u(t)^\t)dt + [ T e(t)<t>(t)dt, a 0 = 1. 

i =0 Jo t '=0 Jo Jo 

(2.3) 

The consequences of (2.3) are: (i) modulating (2.1) with 4>{t) has transferred the 
derivatives of the original data pair {u(t),y(t)]t 6 [0,T]} into derivatives of a chosen 
known function and (ii) the estimation of unknown initial conditions can be 
totally avoided due to (2.2). 

Consider specifically the nth order complex Fourier type of modulating function set 1 : 


= ie- jm ^(e-^ ot — l) n 

(2.4) 

i n-j-m 

(2.5) 

f 

6 

+ 

i 

l-WI 

■i 1 It-i 
II 

(2.6) 

m = 0, 1,2, ..., M 
« = ( ; ) 

(2.7) 


where u> 0 is called the resolving frequency defined as w 0 = ZnjT and T is the time 
interval of the data block. 

Applying the above modulating function set to (2.3) leads to the following regression 
model 

7 ^(m) = 7(m)0 + e„(m), m = 0, 1,2, ..., M (2.8) 

J The equivalence between (2.4) and (2.5)~ (2.7) follows from the binomial expansion. 
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with the regressor row vector: 7 (m) = frftm), ...,7“(m)] in which 


T 

7/ ( m ) = (— i) n_i f 

Jo 


(2.9) 


• = 0,1,2 n; /(() = {»(<) or u(i). ) 


The parameter vector 6 and the model error e n (m) are defined respectively as 


^ — a\ \ 

e = _a L n 

»1 

\ j 

f T 

Cn(m) / 

Jo 

Introducing the following notation: 

Y = (7o y (0),7o V (l),-,7oW) T 

l 7(0) \ 

r = 7(1) 

^ H(M) ) 

c = (e n (0), e n (l), e n (M)) T 


(2.10) 


(2.11) 


( 2 . 12 ) 

(2.13) 

(2.14) 


the relation (2.8) can be rewritten into a vector and complex- valued regression form 


Y = T6 + e 


(2.15) 
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2.2 MFT Under a Stochastic Framework 


2.2.1 Mathematical Antecedent and Modulated White Noise 

Stochastic Calculus in Mean Square Sense 

Due to the involvement of stochastic processes, all the stochastic calculus operators 
including limits, continuity, integration and differentiation in this dissertation are 
presumed to be carried out in the mean square (m.s.) sense (refer to Appendix D for 
details). In light of Appendix D, especially Corollary 1 and Corollary 2, all the above 
established modulating properties and relations hold true as long as { u(t),y(t),e(t)} 
are n th order m.s. differentiable. 

Maximum Likelihood Estimate 

Armed with this prerequisite, some further elaborations on (2.15) can be continued. 
In order to handle the complex regression form more efficiently, we shall examine the 
real regression form first. Later on, in the implementation section of this Chapter, 
the process of converting a complex regression form into a real regression form will 
be scrutinized. For a real regression equation in the form of (2.15), provided the 
equation error e has jointly Gaussian A/^O, W) distribution, the log-likelihood function 
of e can be written as [11] 

C(0,W) = -y ln(2>r) - - - 1(K- TB) t W-'(Y - W) (2.16) 
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If the covariance matrix W is known , then from the necessary condition for maxi- 
mizing the likelihood function,^ = 0, the well known weighted least squares estimate 
of 0 can be written as 

e = (r r w- 1 r)- 1 r r w 1 y (2.17) 

and this estimate will be the maximum likelihood estimate of 0 or the minimax 
entropy estimate [11]. 

Another nice property, which might be used in the system order determination prob- 
lem, can be stated in the following Lemma: 

Lemma 1 (distribution of posteriori cost function) Let the cost function J{0 ) 
be defined as 

J(0) = (y - t t 0) t w~ 1 (y - r T 0 ) 

then its corresponding WLS estimate is 

0 = (r T jy- 1 r)- 1 r T W' r - 1 y. (2.18) 

If the covariance matrix W of the sequence {e„(m)}, m = 0,1,2 is known, 
then for the posteriori cost function J{0), 

J(9) ~ X 2 {M + 1 -n) (2.19) 

where symbol ~ means “obeys” and h is the dimension of the parameter vector 0. 
Proof: Define output error residuals e as 

e = Y -T0. (2.20) 
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Putting (2.18) into the above definition: 


e = [j-r(r T w'~ 1 r) _1 r T w''- 1 ]e (2.21 ) 

Decompose VP -1 = W~ 2 • and apply linear filter W~ 2 to the residual e: 

w~k = w r "i[/-r(r T w _1 r)- 1 r T w r - 1 ]e 

= [i- w- 1 ir{r T w- 1 r)- 1 r T w-'2]-w-k (2.22) 

^ 

p x 

Directly from the above definition of P, we can prove that P is an idempotent matrix 

P 2 = P. 

From Lemma A. 28 [54] we know that the following statements hold for an arbitrary 
idempotent matrix P: 

• All eigenvalues are either zero or one 

• Rank(P)=trP 

Then in our case 

Rank(P) = trP 

= tr[/- ip-2r(r r vp- 1 r)- 1 r T fp-2] 

= M + i -tr[(r T w'" , r)“ 1 r T w , -*w''-*r] 

= M + 1 — h 

where n is the dimension of 6. 

Paraphrasing from Lemma B.8 [54]: 
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Assume X ~ Af(m,W) and set Z = AX 4- B for constant A and B of 
appropriate dimensions. Then Y ~ Af(Am + b,AW A T ). 

Therefore, for the filtered sequence X = W~ $ e, 

X = W~h~ Af(0,I) (2.23) 

J{6) = i T W~H = X T PX. (2.24) 

Paraphrasing from Lemma B.13 [54]: 

Let X be an M-dimensional Gaussian vector, X ~ A/^O, I) and let P be 
an (M|M)-dimensional idempotent matrix of rank M — h. Then X T PX 
is x 2 (Af — n) distributed. 

Directly applying the above Lemma to J(6), we have J(0 ) ~ X 2 (M + 1 - h). The 
proof is complete. 

In most identification problems, the covariance matrix W is not available in advance. 
But if some knowledge or assumptions about the statistics of the regression model 
error t can be imposed beforehand, it is indeed possible, as shown in the future 

sections, to derive an explicit form of W which may or may not depend on the 

parameters 6. 
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Stationary White Gaussian Noise 


Time domain stationary white Gaussian noise is a symbolic process n(t) = W'o(t) 
(where W'o(t) is a Wiener-Levy process) with mean function 

E{n(t)} — 0 

and correlation function 

~ ^1 — ^2) 

where S is the Dirac delta function. From the definition of m.s. derivative listed in 
Appendix D, this white Gaussian noise is not m.s. differentiable at all. However, 
after some dedicated and lengthy mathematical maneuvers, e.g., pp 313~328 in [17], 
a certain justification for its wide usage can be made. 


To conclude our discussion, let us summarize what we have gained by intro- 
ducing the concept of white Gaussian noise. For one thing, it allows us to 
apply the rules of MS calculus even to processes that are not MS differen- 
tiable, and it greatly simplifies calculations involving the Wiener integral. • • 
More importantly, however, the white Gaussian noise represents an idealized 
form of a continuous physical noise just like the Dirac delta function is an 
idealized form of a unit impulse. Thus, whenever we wish to model a physical 
noise that in reality may consist of densely packed narrow impulses of constant 
energy and random polarity, we may reach for a white Gaussian noise as a 
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suitable mathematically tractable idealization. 

— H .J. Larson and B.O. Shubert (pp 327 in [17]) 

Another vivid example of approximating Wo(f) on [0 ,T] with a sequence of pro- 
cesses, W( k \t)]k = 1,2,..., was illustrated in [42], pp 94~97, where it was proved: 
(i) the sample functions of W^ k >(t) are infinitely often differentiable on [0,T] with 
probability one, (ii) the sequence is smooth in the m.s. sense on [0, T], (iii) 

= 0, t € [0, T], (iv) W^(i) is normally distributed, (v) If t\ < t 2 < 
h < ^ 4 , the increments of W^ k \t) on [t x , t 2 ) and [t^t 4 ) are orthogonal for sufficiently 
large k, (vi) As k —+ 00 , W^(t) — ► W{t) in m.s. uniformly in t € [0,T], (vii) 
£{WW(s)W< fc >(<)} -► E{W(s)W{t)} as k -> 00 . Hence W^(t) for a sufficiently 
large k will not only be m.s. infinitely differentiable, but also infinitely close to the 
ideal Wiener- Levy process. In the rest of this thesis, we therefore view that the white 
Gaussian noise is actually defined as rc(<) = k —* 00 . 

Modulated White Gaussian Noise 

From the discussion in 2.1, we can say that modulating a time domain process using 
a set of modulating functions is equivalent to applying a linear transformation to it. 
For Gaussian distributed random processes, one well known fact is that any linear 
operation performed on a Gaussian process results in another Gaussian process [17]. 
In our case, the modulating affect on a white Gaussian noise e(t), t € [0, T], can be 
exhibited in the following Lemma: 
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Lemma 2 (modulated white noise) If e(t) is white Gaussian noise in the con- 
tinuous time domain with E[e(t\) • e(f 2 )] = • hp(ti — t 2 ) an d e„(m) is the result of 

modulating eft ) by the n th order modulating function set <f> m<n {t) defined in (2.6), i.e., 

[T 

Cn(^) — / < Am,n(t)e(^ )dt 

Jo 

then the sequence {e n (m)}, m = 0,1,2 is a stationary discrete Gaussian pro- 
cess; its covariance matrix W e = E(e*e T ) is a banded Toeplitz matrix with bandwidth 

n, each element of which is real and can be expressed as 

(0 1/1 > n 

^e (m , m+0 = E[e n (m)e* n {m + /)] = l ^ ( -i)>.( 2n )! ... < (2.25) 

l X (n-/)!(n+/)! I‘l - " 


Proof: From the definition and Corollary 3 in Appendix D 


E[e n (m)e*(m + /)] = E[J^ e(t x ) • • J q e(t 2 ) • 4>* m +i,n(h)dt2] 

= £ £ E\e{t x ) • e{t 2 )] ■ ■ <f>* m ^ n {t 2 )dUdt 2 . (2.26) 

Noting that E[e(t\) • e(t 2 )] = a 2 • Spifi — t 2 ) and utilizing the sifting property: 

/ 9(t) ■ $D(t ~ r)dt = g(r) 

Jo 

for any continuous function g(t), equation (2.26) can be written as 
E[e n {m)e m n (m + /)] = a 2 f <j> m<n {t)<t>* m + ln {t)dt. 

JO 


(2.27) 


Considering the relations (2.4), (2.5), and (2.6): 


1 ki=0 k 2 = 0 


(2.28) 


Jti =0 fc 2 =0 
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and substituting this into (2.27) 


£MmK(m + ()] = £ t t c tl c k2 f (2.29) 

1 kj=0k 2 =0 J ° 

But only the frequency indices satisfying l + k 2 — k x — Q will contribute to the above 
integral, i.e., 


E[e n (my n (m + l)} = j ^ 


Ht=n CkCk+1 


> n 
< n 


(2.30) 


and from formula 0.156.2 in [13] 


n-l 

E C k C k+l 

k = 0 



(~l) f -(2n)! 
(n — l)\(n + /)! 


(2.31) 


Combining the above two equations, we have 


^(m.m+l) = E {Cn(my n (m + /)] = { ^ (-l)'-(2n)> !!! / ” (2.32) 

l T (n-/)!(n+0! I* I - ” 

Therefore, the covariance matrix is banded with bandwidth n and each element is 
just a real function of l so that the sequence is at least a wide sense stationary (w.s.s) 
process. But for a w.s.s. Gaussian processes, it must be stationary. Equation (2.32) 
also has manifested the fact that the covariance matrix W e is a Toeplitz matrix. The 
proof is complete. 

More generally, for any jointly Gaussian distributed time domain stationary process, 
its modulated sequence would still be a stationary process, except that the bandwidth 
of the covariance matrix is not necessarily equal to the order (n) of the n th order 
modulating function set. 
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2.2.2 WLS Algorithm and Equation Error Signal Model 


Consider the idealized equation error signal model on [0,T]: 

E a„-,y ( °(0 = E 6n-.-u (<) (<) + e(i), «o = 1 (2.33) 

i=0 i=0 

where the equation error e(t) is white Gaussian noise, and the {ad and {b,}, i = 
1,2, are the time-invariant parameters needed to be identified from the input- 

output data pair {u(t), y(t); t € [0, T 1 ] } . Using the same notation, the corresponding 
modulated equation error model can be put into the same regression form: 

y = r0 + c (2.34) 

where the error vector e = [e n (0), ..., e n (M)] T results from modulating the white Gaus- 
sian noise e(t) using the n th order modulating function set { ), m = 0, 1, ..., M}. 
For the error sequence {e n (m)}, based on Lemma 2, its covariance matrix will be 
W e = E(e m t T ) which in this ideal setting is not related to the system parameters or 
input/output data and hence can be written out explicitly in the form of (2.32). With 
this covariance matrix as weighting and the discussion in 2.2.1, a maximum likelihood 
estimate can be obtained and framed into the following WLS/MFT algorithm 2 : 

Algorithm 1 (WLS/MFT Estimate) 

1. Build the weighting matrix W e = E(e*e T ) based on Lemma 2. 

2 For a detailed implementation of this algorithm, please refer to Section 2.3 
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2. Compute the parameter vector from 

Owls = (r T Wf 1 r)- 1 r r W , - 1 y. (2.35) 

Remark: In most cases the variance a 2 of the noise e(t ) is not known even though it 

has no affect on the above algorithm due to cancellation in (2.35). Thus, if we rewrite 
2 

W e = y We, where We is defined only by binomial coefficients and the order n, then 
(2.35) can be adjusted accordingly as 

Owls = (r T W^- 1 r)- 1 r T W^" 1 y. (2.36) 

2.2.3 Measurement Noise Signal Model and AWLS/MFT 
Algorithms 

Measurement Noise Signal Model 

Different from the equation error signal model, the measurement noise signal model 
shown in Figure 2.1 can be characterized as the ideal input /output data pair {u(i),y(t)} 
contaminated with additive white noises v{t ) and n(t). Our goal here is to identify the 
parameters of model H(s) from this contaminated data pair {u(i) +v(t),y(t) + n(tf)}. 
Assume the model H(s) in the time domain is of the differential form 

E a n -,y (0 (<) = E «0 = 1. (2.37) 

»=o t =o 

Then from Figure 2.1 

u(t) = u(t ) - v(t) (2.38) 

y(t) = y(t) — n(t). (2.39) 
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Figure 2.1: Measurement noise signal model. 

Substituting u(t) and y(t) into model (2.37): 

£>„-»“>(<) = £ k.iuM(t) + ±a^(t) - £ in-it-Wft) (2.40) 

t =0 »=0 »=0 »=0 

^ ' 

e(t) 

Unlike the equation error model, e(t) is directly related to the parameters needed to 
be identified. 


Covariance Matrix of Modulated Error Sequence 

Before modulating the above equation, let us observe a general relation for any n th 
order differentiable function f(t): 

r = (-i y ■ r m ■ ^jt)dt 

Jo Jo 

1 rT n + m 

= ^(-ir * / Ei-jkuoyc^me-^dt 

1 J ° k=m 

1 n+m r T 

= j, £ U k ^o) l Ck-m ■ J f(t) e 3 

k=rn v ^ 

F(k) 


i n+m 

(2.41) 

k=m 


Film ) 


= yfHm) 

(2.42) 
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Armed with this relation, the modulated e(t) becomes 


z(m)=f e(t) • (j> mi n(t)dt 

Jo 

=1 a n _,n (,) (i) - 53 

J0 L :=0 i =0 






-1 


er 


(2.43) 


Partition 9 in (2.10) as 




(2.44) 


Then from equation (2.43) and (2.41) 

e(m) = (A f n (m),A^ n _ 1 (m),...,A^ 0 ) ^ j + (K-i (m), ..., V o (m))(-0 6 ) 
= - 53 c k-m (0^o) n ,...,(i^o) 0 ) { 0 JJ o n(t)e~ ikuot dt 


a(k t m, 0 a) 


i n+m j 

+ 5 ; £ c,.„ ( 0 'fc*»)"‘ 1 ,-, 0 'twb)°) (-«») [ v(t)e-i l “*<dt 

k—m ^ 1 1 — — ' v ■ ■ i -— / ® 

2 n+m -J 2 n+m .7- 

= - £ a(k , m, 0 O ) / n(<)e- J ' w <ft + - 53 /?(*, to, 0 6 ) / v{t)e- jk ^ l dt 


(2.45) 


which serves to define the parameter dependent frequency functions a(fc, m, # a ) and 
f3(k,m,9b). If the statistics of n(t) and u(t) axe known, the above relation would be 
the starting point of computing the covariance matrix of the residual error frequency 
sequence {e(m),m = 0, 1,2, In a special case where n(t) and v(t) are mutually 

independent and (approximately) white Gaussian, we have the following Lemma: 
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Lemma 3 Ifn(t) and v(t) are mutually independent white Gaussian noises , then the 
covariance matrix W s of the modulated error sequence {e(m), m = 0,1,2 is 
banded with the order of the modulating functions as bandwidth, and its (m, m + /) 
element denoted as is real and can be expressed as 

«V, m+ ,) = E[e(m)e*(m + l)] 

( 0 |/| > n 

= | +m + l,m,0 a )a*(k 1 + m + l,m + l,9 a ) 1/1 < n 

( E£ 2 =o ^(^2 + m + 1-, Tn, 0(,)/3*(k2 + m + /, m + /, #t,) 

(2.46) 

or in brevity 

W s = g(6) and W,, is real. (2-47) 


Proof: Denoting the first term of E[e(m)e*(m + /)] by W\ s t: 


1 n+m n+m+l fT rT 

w ut =— Y, E (x(k 3 ,m,O a )-<x-(k 4 ,m+l,e a ) ^ y o ^[n(< 1 )-n(i 2 )]e -^ wo<l+ ^“ ot2 ^ 1 ^ 2 


n+m n+m+l 
ks =m ki =m+l 

For the white noises E[n(ti) • n(t 2 )] = — ^ 2 )]? this leads to 


1 n+m n+m + / -y 

= ™E E a(fc 3 ,m,0 a )a*(fc 4 ,m + /,^) / e^-^dt. 


VFi 


1st 


fc 3 =m £4 =7Ti+/ 


Letting k 5 = k$ — m, and k e = k 4 — rn, — /, we have 


fPi 4 t = ^E E a (^ 5 + m,m,0 Q )a*(A: 6 + m + /,m + /,0 a )/ e 

4 l n L a ‘'O 


fc 5 = 0 * 6=0 


Replacing with fcj and noting that the integral is nonzero only for k$ — k & — l = 0, 


w = / ° 2 Kl > " 

Ut \ + m + f,7n,0 a )«*(^i + m -f /,m + /, 6 a ) |/| < n 
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From the independence assumption, the second term of E{t(m)t*(m + /)] can be 
similarly shown to be 


W 2nd 


=1 


> n 


£2 

T 


Edo 2 + m + l,m,0 b )/3 m (k 2 + m + /,m + /, 6 b ) |/| < n 
Based on the definitions of a(k,m,0 a ) and /3(k,m,0 b ) in (2.45), the simple yet im- 
portant facts that 


f a(k, m, 0 a ) ■ a*(fc, m, 6 a ) > 0 
\ f3(k,m,0 b ) ■ 3*(k,fh,0 b ) > 0 


for all m and m (2.48) 


can be easily drawn. Thus, combining and W 2nd with (2.48) yields (2.46) and 
(2.47). This proves Lemma 3. 


AWLS Algorithms without /with Stability Constraint 

Due to the dependence of a( ) and /?(•) in the covariance matrix on the parameters 
desired to be identified, the maximum likelihood estimate discussed in Section 2.2.2 
cannot be implemented directly. From a numerical point of view, we have the follow- 
ing posed problem: 

How to find solutions for 0awls and W a from the implicit nonlinear equa- 
tion set 

0 awls = (r T iF s - 1 r)- 1 r T w J - 1 y (2.49) 

Wa = g(0AWLs)- (2.50) 

To solve it, the following relaxation scheme, which will be referred to here as the 
adaptive weighted least-squares (AWLS) algorithm, can be constructed; it can be 
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viewed as an approximate maximum likelihood estimate and is similar to the algo 
rithm presented in pp 47~50 of [11] in terms of numerical relaxation: 

Algorithm 2 (AWLS/MFT estimate) 3 

1. Choose an error tolerance p for convergence judgement and initialize the step 
index i = 0. 

2. Estimate initial Bawls with identity matrix I or W e in (2.49). 

3. Give an initial estimate ofW s ° from (2.50) where g{0) is defined by (2.45)~(2.47). 

4 . i—i-hl. 

5. Estimate 0 ' AW ls f rom (2.49) with W^” 1 . 

6. Compute Wl from (2.50) with 0\ WLS . 

7. tf Pawls ~ ®'awls II ^ l 1 ) sto Pi otherwise continue. 

8. Go back to step 4. 

Remarks on the above AWLS/MFT Algorithm: 

1. Unlike the residuals in Lemma 2, the sequence {e(m)}, m = 0, 1, ..., M, for the 
residuals in Lemma 3, is no longer stationary. 

2. Because only a biased estimate 9 awls can b e obtained, so likewise the covari- 
ance matrix will be biased. Therefore, the AWLS estimate is no longer the 

3 For a detailed implementation of this algorithm, please refer to Section 2.3 
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exact maximum likelihood estimate. But when the bias is small, it is still close 
to the maximum likelihood estimate. 

A i 

3. The linear filter W7 is no longer a whitening filter as in Lemma 1. But when 
the bias is small, Lemma 1 still holds approximately. 

4. Since implementing the algorithm depends only on the ratio of the pair (cr£, <t£), 
it can still be implemented for unknown noise levels when either a l or al is zero, 
i.e., one is negligible relative to the other, or when a* « cr^. 

5. When the distributions of n(t) and v(t ) are not known beforehand, the above 
algorithm still could be applied and in almost all our cases gives better results 
than the LS/MFT and WLS/MFT algorithms. 

6. Numerical experiments show that if the algorithm converges, it will converge to 
the same value no matter from which initial weighting, i.e., either from LS or 
from WLS, even though the rate of convergence could be different 4 . 

7. The AWLS algorithm is much less sensitive to the chosen modulating bandwidth 
ub, implying that it is a more robust algorithm 5 . 

4 Using the estimate of 6 with the WLS/MFT algorithm to estimate W° is more likely to lead to 
a faster convergence. 

5 This will be further illustrated in 5.4.2. Please refer to Pearson and Lee [39] for guidelines in 
determining the modulating bandwidth wg and the resolving frequency u 0 . 
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8. As for convergence of the above algorithm, the comment from Goodwin and 
Payne (see pp 50 of [11]) regarding their algorithm seems also valid for this 
algorithm: "Unfortunately the authors are not aware of the existence of a global 
convergence proof for the above relaxation algorithm. However, computational 
studies indicate that the algorithm works well in practice." As a matter of fact, 
Algorithm 2 has never failed to converge in all our numerical study examples. 

In most control systems, the estimated model is required to be stable. Inspired by 
the Projection Algorithm on pp 367 of [22], the Algorithm 2 could be easily adjusted 
according to the following form so that the estimated model with stable poles can be 
essentially guaranteed. 

Algorithm 3 (AWLS/MFT estimate with stability constraint) 

1. Choose an error tolerance (i for convergence judgement and initialize the step 
index i — 0. 

2. Estimate the initial 0^ WLS with identity matrix I or W e in (2.49). 

3. Obtain an initial estimate ofW s ° from (2.50) where g(6) is defined by (2.45) (2.47). 
4 ■ i=i+l. 

5. Estimate 0 AWLS from (2.49) with 

6. If 9 awls n °t stable, i.e., the polynomial a n -»s‘ is not Hurwitzian, mirror 

the unstable poles into the left-half-plane and recalculate 0' AWLS ; otherwise skip 
this step. 
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7. Compute W J /rom (2.50) with 0 \wls- 

8. If ||^vvL5 — O'awlsW — 8 stop; otherwise continue. 

9. Go back to step 4. 

Remark: In most cases, this adjustment is not necessary. However, when there exist 
marginally unstable poles, the above algorithm has been proven to be very effective 
(as a typical example, see Section 5.4.2). 

2.3 Implementing WLS/MFT and AWLS/MFT 

As mentioned earlier, both the real and imaginary parts of Definition (2.4) are still 
modulating functions in themselves. Therefore, modulating a differential system will 
provide two sets of algebraic equations corresponding to the modulated real and 
imaginary parts respectively. In this section, we shall devise a joint cost function 
which can utilize both parts and can be easily minimized under the weighted least 
squares framework. But first, some basic relations bridging the real and imaginary 
parts should be explored and disclosed. 

2.3.1 Basic Relation of Covariance Matrices 

Considering the complex modulating function set (2.4) defined as 

= ^e~ jmwot {e~ jwot ~ l) n (2.51) 
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where n is the order of modulating function set and m = 0, 1, 2, A/, it can be 
separated into a real part <f>^ n (t) and imaginary part ^ n (t) : 

K,n(t) = <,(0 +;>!,„(<)■ (2.52) 

Denoting the signal quantities modulated by <j^ n (f) with subscript R and those 
modulated by <Am,nW with subscript /, we have (cf. (2.12) rs j (2.15)) 


Y = Yr+jYj 

(2.53) 

T = Tn + jTj 

(2.54) 

t = tR + jtl- 

(2.55) 

The regression Y = Y0-\-t can also be divided into two parts such that 

(■R = Yr — Tr6 

(2.56) 

t I = Y I - TjO 

(2.57) 

Further, define four different covariance matrices ( W , W R , Wi and Wrj) 6 

with their 

(m, m + l) th element as 

= ■ t*(m + /)] 

(2.58) 

WR(m,m+i) = E[e R (m) ■ e R (m + /)] 

(2.59) 

= E[ei(m) • ej(m + /)] 

(2.60) 

WR/(ni,m+/) = E[eR{m) • e/(m + /)] 

(2.61) 


6 For the equation error signal model, they become ( W e , We B , W ei and W eRI ). Accordingly, for 
the measurement noise signal model, they will be designated by (W,, W*,,, W a , and W Sm ). 
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where m — 0, 1, 2, M, l = 0, ±1,±2, .... Then for the signal models mentioned 
earlier, we have 


Lemma 4 (relations among the different covariance matrices) Under the mea- 
surement noise signal model as shown in Figure 2.1 and the modulating function set 
(2.51), if n(t) and v(t) are mutually independent white Gaussian noises, then the 
following relations hold true: 

_ 1 JO ; m ^ 0 or / ^ 0 

+ { &a 2 n + &b 2 n ■ m = 0 and l = 0 

_1 JO ; m ^ 0 or / / 0 

U,S/(m ' m+ 0 “ l a 2 n + £bl ; m = 0 and 1 = 0 

where w 4(m m+|) is defined in Lemma 3, eqn. (2.46), and 


(2.62) 

(2.63) 


Ws RI(m,m+r) O’ 


Proof: From definition (2.55) 

eft(m) = ^[e(m) + c*(m)] 


therefore 


(2.64) 


(2.65) 


^R( m,m + f ) = E[e R (m) • e R (m + /)] 

= {[ e ( m ) + e *( m )] ‘ [ e ( m + 0 + e*(m + /)]} 

= ^ E {e(m)e(m + /) + e(m)e*(m + /) + e*(m)e(m + /) + e*(m)e*(m + /)} 

(2.66) 
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Here we need first to prove the fact that terms like £[c(m)e(m + /)] and E[e*(m)e*(m-\- 
/)] will contribute to equation (2.66) only when m = 0 and / = 0 due to the orthogo- 
nality of {e~i ku>ot , k = 0, 1, 2, ...}. 

Using the condition of mutual independence, equation (2.45), and Corollary 3 in 
Appendix D, 


2 n+m n +m- f / vj 

E[e(m)e(m+l)] = ^ 5Za(A? 1 ,m,^a(A:2,m+/,fl 0 )/ E[n(t 1 )n(t 2 )]e~^ klWot+kiW0 ^dtidt 2 + 

ki=mk2=m+l ® ® 


T tT 


IIi 


2 n+m n+ro+/ -j- 

7 ~ E E^(*i^A)/?(*2,m+/A)/ /£;[u(f 1 )^(^)]e- j[<:iW0i+i2W0,] df 1 ^ 

= m ^2 — m +/ ^ ^ 


T rT 


n 2 


(2.67) 


Note the fact that 


E[n(ti) ■ n(t 2 )} = <t^d(<i - f 2 )- 


Then the first term of equation (2.67) can be reduced to 

r 2 n+m n+m+J -j 1 

m, 0 a )a(A: 2 , m + l,0 a ) j 

k\ =m k2 = m+l 

In the above, m > 0 and m + / > 0, so Aq > 0 and k 2 > 0. But due to 


2 n+m n+m+/ .j 1 

n i = ^E E cc{k u m,0 a )cc{k 2 ,m + l,0 a ) e ~^ +k ^ 'dt (2.68) 

-I 1 1 1 / ‘'O 


J T e~ j{kl+k3)wot dt = | 


we have 


n 1 = 


0 


0 ? fci -h /^2 ^ 0 

T j Aq H~ ^2 = 0 1 


; fci + /?2 7 ^ 0 


^ ?p-a 2 (0, 0, 6 a ) ; ^ = 0 

Condition ^ -f k 2 = 0 is equivalent to m = 0 and l = 0. Further from the definition 

a(k, m, 6 a ) = c k - m (O'Am,)", (jku>o)°) ^ ^ (2-71) 


(2.69) 


(2.70) 
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we have a(O,O,0 a ) = —a n . Therefore, 


jj_fO ; m^Oor/^O 

1 \ ; m = 0 and / = 0 

Following a similar argument, we have 

j-j / 0 ;ra^0or//0 

i T^n ) rn — 0 an< l / = 0 

as well. Substituting III and n 2 into (2.67) 

£[e(m)e(m + /)] = ( ^ d ’ 

l T U n ' T U n i 

is proved. This also automatically implies 


m ^ 0 or / / 0 
m = 0 and l = 0 


(2.72) 


(2.73) 


(2.74) 


£[e*(m)e*(m + /)] = {E[t(m)e(m + /)]}* = j 




£i. b 2 

T u n 


; m / 0 or I ^ 0 
; m = 0 and / = 0 

(2.75) 


Using (2.74), (2.75) and Lemma 3, equation (2.66) can be rewritten as 

W *R(m.m + l) = 2 ' {■^[ C ( m ) e *( m + 0] + - E [ C ( m ) C ( ? ^ + 0]} 

_ 1_ f 0 ;m^0or/^0 

2 | -f- ; m = 0 and / — 0 


(2.76) 


If it is further noted that 


£ /( m ) = i[e(m) - c *( m )] 


(2.77) 


then 


tD */ ( m.m + o = E[ti{m)ti{m + /)] 
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— ■ E{[e(m) - c*(m)][c(m + /) - e*(m + /)]} 

— E{e(m)e(m + l) — e*(m)e(m + /) — e(m)e*(m + 1) + e*(m)e*(m + /)} 


Similarly considering the results of (2.74), (2.75) and Lemma 3, the above relation is 
reduced to 


t) 


- • {£[e(m)e*(m + /)] - E[e(m)e(m + /)]} 

1 f 0 jm^Oorl/O 

- W a 2 + |i 6 2 . m = o and / = 0 


(2.79) 


Hence, relations (2.62) and (2.63) are established. Finally, using identities (2.65), 
(2.77), (2.74), (2.75) and Lemma 3: 


5 JU(m,m+l) 


= E[tji(m)ei(m + /)] 


= — • E{[e(m) + e*(m)] • [e(m + /) - e*(m + /)]} 

V 

= l - ■ Irn{w S(m m+l) } + jj- E[e(m)e(m + /)]- — E[t{m)e(m + /)] 

= 0 (2.80) 


Therefore, the proof of this Lemma is completed. In the form of a matrix, the above 
Lemma implies: 




(2.81) 


W.* = 2 W + Mf). 


(2.82) 


= ^{w. - wJ), 


(2.83) 
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where the column vector r) is defined as 



As a matter of fact, Lemma 2 and the equation error signal model could be just 
deemed as a special case of Lemma 3 and Lemma 4, and correspondingly, t] e should 


be modified as: 



(2.85) 


and 


Wc R J ~ 0, 

(2.86) 

), 

(2.87) 

w„ = \{W t - V , V J). 

(2.88) 


Hence, the following derivations about implementation schemes will be applicable to 
both stochastic signal models, and subscripts “e” and “s” will be dropped for the 
general discussion. 


2.3.2 Joint Cost Function 

Splitting the modulated quantities into real and imaginary parts is equivalent to 
using both real and imaginary parts of the complex modulating function set (2.4) to 
modulate the original differential equation model separately. Although it would be 
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quite awkward, both real and imaginary parts could work independently to estimate 
the very same parameter vector 6 using the weighted least square algorithms with 
different weights Wr and Wi through minimizing two different cost functions. In 
order to avoid the potential agony induced by this separation, the possibility of using 
one joint cost function binding these two parts to produce just one estimate, instead 
of two, should be explored. Under the assumption that additive noises n(t) and v(t) 
are mutually independent white Gaussian noise, the error sequences tR and tj in 
(2.56) and (2.57) are uncorrelated as shown in Lemma 4, and then the joint likelihood 
function is as follows: 


p{tR,£i | 9) = p(cr I 0)-p{ti \ 9) 

1 


(2t)"+‘|H / r|!|H / ;| 5 ' exp {~2' / ^} <2 ' 89! 


where 


J{9) = elW^tn + eJWf 1 ^ 

= (y r - r^) 7 ^ 1 ^ - r*0) + (Yi - Tjefwr'iYj - r,<?) 
> o. 


The log-likelihood function C{9) can be used and 


£(0) = lnp(e/?, ej\9) 

= -(M + 1) In 2tt - |[ln \W R \ + In | W,|] - l -J{9). (2.90) 
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If the covariance matrices Wr and Wi are known, maximizing C{9) for the maximum 
likelihood estimate is equivalent to minimizing the quadratic form of the function 
J(0). For this reason, 

J(9) = (Y r - V R 0) T Wn\Y R - r r B) + (Yj - TrffWf'iYj - Tj8) (2.91) 


has been selected as the joint cost function for the rest of our studies on SISO systems. 
The combined or joint estimate desired should be 

8 = arg min J(8). (2.92) 

From the necessary condition of minimization: |^ = 0, we have 

% = r T R w i ;'(Y R -r R e) + r';wf'(Y l -r I 9) 

= (r*tVUi + Tjwr'Y,) - (r5w£‘r„ + rjw7 1 r,)0 

= 0 (2.93) 


which implies 


i = {r^T* + rfiypr,}-' • {1^% + tfwr'Y,}. (2.94) 

In order to have a WLS form like (2.35), further introduce the following combined 
notations: 


r c = 

(?:) 

(2.95) 

Y c = 

(K) 

(2.96) 

W c = 

( W R 0 \ 

0 Wj ) ■ 

(2.97) 
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If we denote 



w *- 1 0 \ 

o W7 1 j 


then it is straightforward to show that 




= (rJtyjTcj-Tj^'Kc 


-1 rT w-li 


(2.98) 


This result shows that by minimizing the joint cost function (2.91) we still can have a 
combined WLS estimate 0 through the combined regressor, regressand and weighting. 
This also has provided an efficient way to utilize the information carried in both the 
real and imaginary parts of modulated quantities. 

Further, we have the following observations about implementing the above scheme: 
1. For the regular least squares algorithm, Wc = /, then 

ks = (r£r c )-'rjy c . 


2. For the equation error signal model and measurement noise signal model of 
Figure 2.1, two matrix inverses Wr 1 and VFj -1 are needed for W ^ 1 . However, 
utilizing the matrix inversion lemma [2]: 

(A + BCD)- 1 = A' 1 - A^BiC' 1 + DA~ 1 B)- 1 DA~\ (2.99) 

and letting A = VF, B = r), C = I and D = »/ T , the special form (2.82) of Wr 
becomes 

Wr 1 = 2 (W A t/t/ 7 )- 1 = 2(VF -1 - W~ l ri(l + t fW-'xiY'rFW- 1 ). (2.100) 
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j 2 2 

Denoting S 3 = y for the measurement noise signal model (or 6 e = 

yfif for the equation error signal model) and A = t?t ] T W~ 1 (with a similar 
designation A s and A e ), and partitioning 


V = 


6 

OmxI 


we have 


W~ l = ( WlU W [12 \ 

' V ) ’ 


( 2 . 101 ) 


t fW-'-q = S 2 wm 


A = rjr] T W- 1 

( S 2 w m 6 2 W m \ 

\ 0\lx 1 OmxM j 


(2.102) 


(2.103) 


It is straightforward to show 


W* 1 


2W~ 1 (I 


A 

1 + 6 2 wiu 


)• 


(2.104) 


Similarly, we have 

= m " {l + rdb^)- ( 2105 > 

Equations (2.104) and (2.105) indicate that only one matrix inverse VK -1 is 


needed at each iteration. 


3. One very important remark that should be reiterated here is the fact that if 
W R and Wj are not known beforehand, while they may be explicitly expressed 
as a function of the parameter 0 such as shown in Lemma 3 and Lemma 4 
for the measurement noise signal model, the AWLS/MFT estimate stated in 
Algorithm 2 does not lead to the exact maximum likelihood estimate. But 
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when the estimated 9 through minimizing J(0 ) is not far away from the true 9 , 
the estimated W should also be very close to the true one as well. This is why the 
claim that AWLS/MFT is just an approximated maximum likelihood estimate 
has been declared. In order to have a true maximum likelihood estimate, the 
maximization should be applied directly to the log-likelihood function (2.90). 
The terms like In |W«| and In \Wj\ could truely make the computation become 
formidable. In this aspect, J{9) is much more attractive. 


4. For computational simplicity, one may simply use a single (approximate) real 
W as the weighting for both W R and Wj in each iteration and neglect the true 
updating forms (2.104) and (2.105) caused by the tiny difference between W R 
and Wi. As one revealing example to acertain the cost of this simplicity, let us 
identify the following second order system: 


H(s) 


8 

s 2 + 4s -f 10 


(2.106) 


The AWLS/MFT Algorithm 2 will be utilized with and without forms (2.104) 
and (2.105) using 100 Monte Carlo simulation runs for each case at each of 
several additive noise levels. In order to see the relative difference one to another 
on both the mean and standard deviation, define a percent error measure by: 


A = ,f a|12 • 100%. (2.107) 

Here £ e corresponds to either the mean or the standard deviation values for 
each parameter obtained with the exact weightings as defined in (2.104) and 
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(2.105), and £ a corresponds to either the mean or the standard deviation values 
for each parameter obtained without the exact weightings as defined in (2.104) 
and (2.105). The results are summarized in Table 2.1, from which we can see 


true para 

1 __J 

r 4 

10 



mean 

std 

mean 

std 

mean 

std 

NSR 

EXACT 

7.9920 

0.08105 

4.0025 

0.04378 

9.9970 

0.07880 

5% 

APPROX 

7.9862 

0.08700 

3.9997 

0.04579 

9.9930 

0.08265 

EXACT 

7.9661 

0.15862 

3.9795 

0.09232^ 

9.9710 

0.1424 

10% 

APPROX 

7.9621 

0.16976 

3.9777 

0.09795 

9.9688 

0.1513 

EXACT 

7.9033 

0.30242 

3.9500 

0.17236 

9.9160 

0.29532 

20% 

APPROX 

7.9032 

0.33463 

3.9503 

0.18710 

9.9164 

0.31104 

EXACT 

7.5696 

0.61485 

3.7424 

0.33769 

9.6179 H 

0.54303 

40% 

APPROX 

7.5373 

0.61791 

3.7259 

0.34017 

9.5889 

0.55101 

EXACT 

6.6254 

0.95782 

3.1916 

0.54162 

8.7639 

0.78983 

80% 

APPROX 

6.5373 

1.00141 

3.1437 

0.55088 

8.6856 

0.79628 

A 

0.55% 

4.68% 

0.60% 

2.78% 

0.39% 

2.08% 



100 Monte Carlo runs: 

APPROX means using Ws as both Ws R and Wsj without using (2.104) and (2.105). 
EXACT means using exact the Wsr and Wsj as defined in (2.104) and (2.105). 


Table 2.1: Comparison between using exact and approximate weighting matrices 


that the increased accuracy in using (2.104) and (2.105) has only a slight edge, 
i.e., 0.6% in mean and 4.68% in std, over the case without using (2.104) and 
(2.105). As for the speed of convergence and computational time concerns, they 
do not exhibit any difference. In the rest of this thesis, the results in applying 
AWLS/MFT are obtained without using (2.104) and (2.105). 
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2.3.3 A Simple Recursive Algorithm for Matrix Inversion 


We have seen that one of the major computational burdens is the inversion of an 
(M + 1) x (M + 1) weighting matrix needed in implementing the WLS/MFT and 
AWLS/MFT algorithms. Numerical experiments show that when the order of the 
modulating function set or the model goes higher than 10, the MATLAB’s matrix 
inverse routine which uses “matrix division” or singular value decomposition would 
fail to provide usable answers. Part of the reason is that those routines are not 
specifically written to deal with banded symmetric positive definite matrices like our 
weightings. The round-off errors could accumulate very fast or the matrix could be 
badly scaled, especially when M gets large. With these particular sparse matrix 
structures in mind, we hope that we can contrive some algorithm which eventually 
avoids direct matrix inversion and also can utilize the sparse structure of the weighting 
matrix to improve the numerical accuracy and efficiency. 

If the upper-left (k + 1) x (k + 1) sub-matrix of the weighting W is denoted by Wk+i, 
we need first to answer the following question: 

Provided that VF fc -1 is known, is it possible to compute from VF fc -1 

without employing a multidimensional matrix inversion ? 

Partition the (& + l)x(A; + l) sub- weighting matrix VFj t+i 

w »' - ( W bi ll ) ( 210S > 

where Bk is k x 1 column vector and a* is the ( k + l)-th diagonal element of the 
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weighting matrix W . Due to the symmetric positive definite property, the inverse of 
W k + 1 exists and is also symmetric positive definite. If we denote 


W -1 ( Vk c k \ 

Wm = l Cl * j 


(2.109) 


where V k is a k x k matrix, C k is k x 1 column vector and d k is the (k + l)-th diagonal 


element of W 1 , it follows from the condition 


Wib+i ' w k+i 


( w k B k \ ( V k C k \ 

V a k j \ C k d k J 

( W k V k + B k Cl W k C k + B k d k \ 

V B T k V k + a k C T k B T k C k + a k d k j 



that we have four equations: 


Bl C k + a k d k 

= 1 

(2.110) 

W k C k + B k d k 

= 0 

(2.111) 

B T k V k + a k C T k 

= 0 

(2.112) 

W k V k + B k C T k 

= h. 

(2.113) 


We wish to solve for the three unknowns (V k ,C k ,d k ) with ( W k , W k x ,B k , a k ) as knowns. 
From equation (2.113) and (2.110) we have 


= W k ~' ■ (I k - B k C T k ) 
, 1 - B\C k 

d k = . 

a k 

Substituting (2.115) into (2.111): 

W k C k + B k l ~ B *° k = 0 

dk 


(2.114) 

(2.115) 
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(Wk- 


B k Bfi _ Bk 

ax a k 


So 


c k = -[W k - ?!%-]-' • ^ 


ak 


ak 


(2.116) 


This expression for Ck still requires computing a k x k matrix inverse [Wk 0fc *'] 
which is not desired. But from the well known Matrix Inversion Lemma, 

BtBL., B k 


C k = -[W k --^] 


ak 


ak 


= -[W * -1 - W^B k iBjW^B k - a k )-'B T k W; l \ 


-11 Bk 


a k 


— — [Ik + 


W k 'B k Bj W;'B k 
a„-BlW k -'B t ‘ a k 


W k 'B„„ BlW^Bj, 


ak 

W k -'B k 


: ri _ tUi-LLh — \ 

' L BlW k x B k - a k 


B T k W k l B k -a k 


Putting this back into equation (2.115) and (2.114) we obtain 

1 


dk = 


' B T k W k l B k -a k 


V - W-H, _ MQg 1 

Vt ~ Wt |7 ‘ B k W k l B k — at 
Further, by defining a column vector A k hy A k = W k 1 Bk and a scalar r k = B k W k 1 B k , 
and combining (2.117) with (2.119) and (2.118), we have the following algorithm for 
any positive definite matrix W. 


(2.117) 


(2.118) 

(2.119) 


Algorithm 4 (recursive matrix inversion for general weighting) For a known 
M xM positive definite matrix W with main diagonal elements {a^}; k = 0, 1, ..., M — 
1, its inverse can be computed in the following way: 
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1. Initialize index k = 1 and W, 1 = — . 

1 O0 

2. Increment k=k+l and check if k > M or not ? Yes: stop; No: continue. 

3. Obtain Bk and ak directly from partitioning Wk as in the right side of (2.108). 
4- Compute 


5. Form 


Ajt 

Tk 

C k 

Vk 

dk 


WF 1 ■ B k 
B T k ■ A * 

Afc 

Tk ~ Ok 

w -i _ AjAL 

n ~ a k 


1~k &k 




( v k c k 

\ Ck dk 


) 


6. Go to step 2. 


Therefore, the question of recursive updating W k ^ from W k 1 has been answered. 
Following up, we will further take the sparse structure of the weighting matrices into 
account, so that the above algorithm can be made more efficient to compute A k an d 
Tk at each recursion. If the order of the modulating function set is n, from Lemma 
2 and Lemma 3 the bandwidth of the covariance would be n. Therefore, except for 
k < n, the column k x 1 matrix can be partitioned as 


B k 


I Q{k-n)x 1 

{ 


) 


(2.120) 
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where £)( fc - n ) xl is a (k — n) x 1 zero column vector and B k xl is a n x 1 matrix. 


Correspondingly 


/ ry(k—n)x(k—n) rj{k~n)xn \ 

W-l _ I "mi U k,12 \ 

‘ ( (4^ ,X *) T «ijS )' 

(2.121) 

Then it is straightforward to obtain that 


a _ l Rku )xn ■ sr 1 \ 

^k 1 T^nxn rynxl I 

V K k , 22 * &k ) 

(2.122) 

r„ = ffii'fntSBt* 1 

(2.123) 


Armed with these two definitions and Algorithm 4 we are ready to introduce the 
recursive algorithm for a banded symmetric positive definite matrix. 

Algorithm 5 (recursive matrix inversion for banded weighting) For a known 
M x M banded symmetric positive definite matrix W with bandwidth n and main di- 
agonal elements {a k };k - 0, 1, ..., M- 1 , its inverse can be computed in the following 
way: 

1. Initialize index k = 1 and Wf 1 = — . 

4 Cl0 

2. Increment k=k+l and check if k > M or not 9 Yes: stop; No: continue. 


3. Obtain B k and a k directly from partitioning W k as in the right side of (2.108) 
and (2.120). 


4. Partition W k and B k as (2.121) and (2.120) and compute 


A* 


D (k-n)xn nnxl 
n k, 12 _ ' D k 

nnxn pnx 1 
* D k 
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5. Form 


Tk 


Ck 


Vk 


dk 




nXn pnX 1 
k , 22 D k 


A k 

Tk ~ dk 



Tk ~ dk 


1 


Tk ~ dk 





e* 

4 


) 


Go to step 2. 


Remarks on the above algorithm: 

1. As only the last n columns of W are involved, the maximum inner product 
dimension is n instead of k. When k » n, this is very helpful for the depression 
of accumulation errors. Computationally, kn flops for A* and n 2 flops for Tk are 
required to update at each recursion. 


2. The most computationally-demanding term in updating Vk is A*Aj ; however, it 
only involves the product operations among the elements of the column vector 
A*. Hence, it does not contribute to accumulation errors at each recursion. 

3. The total flops required is of order 0(M 3 ), which is the same as LU decompo- 
sition and Gauss-Jordan elimination methods. 


4. When W is Toeplitz, i.e., W = W e , we have two choices: 
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(a) Still use Algorithm 5. 


(b) Employ Trench’s algorithm (see pp 132 of [9]) which only needs 0(M 2 ) 
flops. The Trench Algorithm requires the column vector Cm and d,M to 
be obtained first from Durbin’s algorithm which is of 0(M 2 ) flops. Note 
in Algorithm 5 that if we drop the V k updating, we also can obtain Cm 
and d,M in 0(M 2 ) flops which is as efficient as Durbin’s algorithm. Then 
at least we can use the above algorithm to first obtain Cm and 

5. The bottom line is that numerically the above algorithm is much more robust. 
It has successfully inverted a M = 1024, n = 12, Toeplitz matrix while the 
routines in MATLAB failed. For AWLS/MFT, it is as efficient as any other 
inversion algorithm. 

6. One other by-product of the above algorithm is that it facilitates writing a 
recursive weighted least squares algorithm, which might not be necessary in 
MFT, but it may be of value to other sequentially correlated data analyses. 


2.4 Comparing LS, WLS and AWLS with PEM 

The second order system: y(t)-\-3y(t)-\-8y(t) = bu(t), where 0 < t < T and T = lOsec, 
was used to evaluate and compare the performance of the LS/MFT, WLS/MFT 
and AWLS/MFT algorithms and to compare with a commercially available PEM 
(prediction error method) algorithm [20] in MATLAB written by L. Ljung [21]. Two 
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hundred Monte Carlo runs were made at each of several noise-to-signal ratios for 
additive white output noise corrupted data. The input signal was u(t) = sin(t 2 /5), 
t € [0, 10] secs for each run and the sampling rate is fixed as 25.6Hz. The output y(t) 
is a combination of the simulated output using LSIM() of MATLAB and the white 
Gaussian random noise sequence generated by RANDNQ, i.e., 

y{t) = LSIM (A,B,C, D,u,t,X0) + n(t) 

where [A,B,C, D] = TF2SS (5, [1 , 3, 5] ) , X0 is the initial condition, and n(t) = k * 
RANDN(256, 1) is the additive noise with the scale factor k determining the noise 
level. In order to have a fair and accurate comparison, every noisy input/output 
realization pair has been forced to run through all four algorithms in each Monte- 
Carlo trial. The noise-to-signal ratio (NSR) is defined as 

NSR = ■ 100% (2.124) 

Ily(0li2 ' 

which characterizes the percent additive noise on the output. As for a true parameter 

0o and its estimate g (with standard deviation <x), a normalized bias and standard 

deviation are formed as 


Normalized Bias = |- — — | - 100% 

0o 

Normalized STD = | — | • 100% 

0 


(2.125) 

(2.126) 


These will be used to measure the accuracy of the different algorithms. For the above 
specific system, its step response will take about 4 seconds to reach steady state and, 
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therefore, for the total 10 seconds of data, the initial condition AO and the input 
could both play an important role. In order to have a better picture of the impact 
of initial conditions on the estimation, two cases with and without randomized AO 
have been carried out in the following simulation studies. 

2.4.1 With X0 Fixed as (0,0)' 

In this case we assume that the initial conditions are always known as .AO = (0,0) , 
so as not to treat the A0 as an unknown in PEM. In each Monte Carlo run of the 
PEM, the initial guess of the parameters is set favorably to the true values as well 
as giving it the true value A0 = (0,0)'. Under this relatively ideal setting for PEM, 
the simulation results are summarized in the Table 2.2 ~ 2.5 and Figure 2.2 from 
which we have the following observations: 

1. Although the PEM has a smaller variance than LS/MFT at most noise levels, 

especially in the lower noise level cases 1 * * * * * 7 , PEM does have greater variance than 

both WLS/MFT and AWLS/MFT algorithms at all the additive noise levels. 

The variances in WLS/MFT or AWLS/MFT have been reduced to about one- 

third the variance of LS/MFT. Between the standard deviations of WLS/MFT 
and AWLS/MFT, the latter has a slight edge over the former only at very large 
additive noise levels. 

7 Fullerton, A. Jr. revealed this fact from his early simulation studies [8]. Our craving of further 

curbing this quantity triggered our studies on the WLS/MFT and AWLS/MFT algorithms. 
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200 Monte Carlo runs for PEM without estimating X(0), Fb=0.4Hz, fixed initial X(0)=(0,0)' 


true parameters 

8 

3 

5 

lly(t)llz=6.472 

mean 

8.509 

3.151 

5.320 

lle(t)llr=0.822 

variance 

0.0171 

0.0069 

0.0153 

NSR=12.7% 

mean 

8.504 

3.155 

5.309 

116(0112=1.664 

variance 

0.0775 

0.0294 

0.0625 

NSR=25.7% 

mean 

8.509 

3.146 

5.328 

Ile(t)ll2=3.330 

variance 

0.2695 

0.1119 

0.2435 

NSR=51.5% 

mean 

8.609 

3.201 

5.417 

1^(0112=4.938 

variance 

0.6306 

0.2530 

03162 

NSR=76.3% 

mean 

8399 

3.181 

5.412 

lle(t)ll 2 =6.660 

variance 

1.3089 

0.5470 

1.1467 

NSR= 102.9% 


Table 2.2: PEM with fixed .Y(O) = (0,0)' and initial guess at true values. 


200 Monte Carlo runs for LS/MFT algorithm, Fb=0.4Hz, fixed initial X(0)=(0,0)' 


true parameters 

8 

3 

5 

Hy(t)H 2=6.472 

mean 

7.991 

2.988 

4.983 

Ile(t)ll2=0.822 

NSR=12.7% 

variance 

0.0294 

0.0126 

0.0290 

mean 

7.961 

2.960 

4.947 

Ile(t)ll2=1.664 

NSR=25.7% 

variance 

0.1350 

0.0557 

0.1283 

mean 

7.886 

2.872 

4.856 

Ile(t)ll2=3.330 

NSR=51.5% 

variance 

03144 

0.1907 

0.4499 

mean 

7.648 

2.698 

4.614 

lleCt)ll2=4.938 

NSR=76.3% 

variance 

0.8119 

0.3434 

0.7714 

mean 

7.077 

2.347 

4.111 

HeCOII 2=6 .660 
NSR= 102.9% 

variance 

1.1768 

0.4449 

1.2680 


Table 2.3: LS/MFT with fixed X(0) = (0,0)' and F\, = 0.4Hz. 
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| 200 Monte Carlo runs for WLS/MFT algorithm, Fb=0.4Hz, fixed initial X(0)=(0,0)' || 

true parameters 

8 

3 

5 

lly(t)ll 2 =6.472 

mean 

7.995 

2.989 

4.985 

116 ( 0112 = 0.822 

variance 

0.0135 

0.0049 | 

0.0103 

NSR=12.7% 

mean 

7.977 

2.967 

4.955 

Ile(t)ll 2 = 1.664 

variance 

0.0580 

0.0228 

0.0443 1 

NSR=25.7% 

mean 

7.859 

2.850 

4.808 

Ile(t)ll2=3.330 

variance 

0.2129 

0.0877 

0.1649 

NSR=51.5% 

mean 

7.696 

2.665 

4.582 

Ile(t)ll2=4.938 

variance 

0.4740 

0.1731 

0.3920 

NSR=76.3% 

mean 

7.306 

2.405 

4.253 

lle(t)ll 2 = 6.660 

variance 

0.4677 

0.1779 

0.3879 

NSR= 102.9% I 


Table 2.4: WLS/MFT with fixed X{0) = (0,0)' and F b = 0.4Hz. 


| 200 Monte Carlo runs for AWLS/MFT algorithm, Fb=0.4Hz, fixed initial X(0)=(0,0) |j 

true parameters 

8 

3 

5 

Ily(t)ll2=6.472 

mean 

7.994 

2.990 

4.986 

116 ( 1 ) 112 = 0.822 

variance 

0.0125 

0.0045 

0.0098 

NSR=12.7% 

mean 

7.974 

2.972 

4.959 

Ile(t)ll 2 = 1.664 

variance 

0.0542 

0.0212 

0.0428 

NSR=25.7% 

mean 

7.863 

2.870 

4.833 

Ile(t)ll2=3.330 

variance 

0.2052 

0.0865 

0.1620 

NSR=51.5% 

mean 

7.703 

2.718 

4.630 

Ile(t)ll2=4.938 

variance 

0.4291 

0.1678 

0.3523 

NSR=76.3% 

mean 

7.338 

2.477 

4.315 

Ile(t)ll2=6.660 

variance 

0.4079 

0.1610 

0.3552 

NSR= 102.9% 


Table 2.5: AWLS/MFT with fixed X(0) = (0,0)' and F b = 0.4Hz. 
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2. As for the bias concern, LS/MFT, WLS/MFT and AWLS/MFT have obviously 
smaller bias than PEM in the lower noise cases, i.e., NSR < 56%. Among 
LS/MFT, WLS/MFT and AWLS/MFT, the bias of AWLS/MFT is smaller 
than that of LS/MFT at the high noise end. 

Overall, under this extremely ideal setup favorable to the PEM, WLS/MFT and 
AWLS/MFT have performed significantly better than PEM within the moderate 
additive noise range in the sense of variance and bias. 

2.4.2 With Randomized XO 

In this case, each realization of an input /output pair is implemented with random- 
ized initial conditions XO. We assume that the initial conditions of the system are 
unknown beforehand, so that XO has to be estimated in the PEM algorithm while 
there is no difference to the MFT algorithms. For each Monte Carlo run of PEM, 
the initial guess of parameters was still set to the true values, but with a randomized 
initial guess of XO which has to be estimated by the PEM in the end. Under this 
relatively thorny condition for PEM, the simulation results are summarized in the 
Table 2.6 ~ 2.9 and Figure 2.3 from which we can make the following remarks: 

1. PEM has not only failed to achieve a smaller variance than LS/MFT, but also 
exhibits unreliability with its frantic-looking mean values. Meanwhile more im- 
portantly, there were almost no noticeable effects on the results of the LS/MFT, 
WLS/MFT and AWLS/MFT algorithms at all noise levels. This kind of robust- 


53 



200 Monte Carlo runs for PEM while estimating X(0) with randomized initial X(0) jj 

true parameters 

8 

3 

5 

lly(t)llj=6.472 

mean 

9.273 

3.561 

5.901 

lie(t)llj=0.822 

variance 

0.6927 

0.3744 

0.7540 

NSR=12.7% 

mean 

10.144 

4.291 ' 

6.904 

116(0112=1.664 

variance 

2.9970 

1.1953 

2.4566 

NSR=25.7% 

mean 

9.102 

3.301 

5.610 

Ile(t)ll2=3.330 

variance 

3.4521 

15663 

2.9923 

NSR=51.5% 

mean 

8.927 

3.201 

5.508 

Ile(t)ll2=4.938 

variance 

4.4852 

1.1124 

2.8240 

NSR=76.3% 

mean 

8.608 

3.145 

5.339 

Ile(t)ll2=6.660 

variance 

5.2385 

0.8092 

2.6997 

NSR= 102.9% I 


Table 2.6: PEM with randomized X(0) and initial parameter guess at true values. 


200 Monte Carlo runs for LS/MFT algorithm with Fb=0.4Hz and rando miz ed initial X(0) 


true parameters 

8 

3 

5 

|| y (t)fl 2 =6.472 

mean 

8.008 

3.000 

5.002 

Ile(t)ll2=0.822 

variance 

0.0225 

0.0081 

0.0201 

NSR=12.7% 

mean 

7.948 

2.959 

4.942 

Ile(t)ll 2 = 1.664 

variance 

0.0923 

0.0334 

0.0932 

NSR=25.7% 

mean 

7.790 

2.865 

4.787 

Ile(t)ll2=3.330 

variance 

0.3036 

0.1131 

0.2717 

NSR=51.5% 

mean 

7.488 

2.708 

4577 

lle(t)Jl2=4.938 

variance 

05972 

0.2411 

0.6038 

NSR=76.3% 

mean 

7241 

2.469 

4243 

Ile(t)ll2=6.660 

variance 

1.0623 

0.4029 

1.02% 

NSR= 102.9% 1 


Table 2.7: LS/MFT with randomized X(0)' and Fb = 0.4Hz. 
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200 Monte Carlo runs for WLS/MFT algorithm with Fb=0.4Hz and randomized init ial X(0) 

5 


true parameters 


mean 


variance 


mean 


variance 


mean 


variance 


mean 


variance 


mean 


variance 


8 


7.981 


0.0125 


7.933 


0.0392 


7.901 


0.1304 


7.683 


0.3013 


7.529 


0.5969 


3.005 


0.0036 


2.977 


0.0125 


2.965 


0.0471 


2.834 


0.0936 


2.664 


0.1792 


4.9% 


0 . 00 % 


4.950 


0.0302 


4.919 


0.09% 


4.739 


0.2541 


4554 


0.4907 


lly(t)ll 2 =6.472 

116 ( 0112 = 0.822 

NSR=12.7% 


Ile(t)ll 2 = 1.664 
NSR=25.7% 


Iie(t)ll2=3.330 

NSR=51.5% 


Ile(t)ll2=4.938 

NSR=76.3% 


He(t)ll 2 = 6.660 
NSR= 102.9% 


Table 2.8: WLS/MFT with randomized X(0)' and Fb = 0.4Hz. 


200 Monte Carlo runs for A WLS/MFT algorithm with Fb=0.4Hz and randomized initial X(Q) 


true parameters 


mean 


variance 


mean 


variance 


mean 


variance 


mean 


variance 

mean 


variance 


8 


7.989 


0.0111 


7.941 


0.0347 


7.890 


0.1229 


7.668 


0.2725 


7.497 


0.5456 


3.006 


0.0034 


2.981 


0.0109 


2.960 


0,0442 

2.850 


0.0857 


2.673 


0.1621 


4.999 


0.0093 

4.955 


0,0273 

4.912 


0.0929 

4.748 


0.2313 

4.538 


0.4479 


|| y (t)ll 2 =6.472 



116 ( 0112 = 0.822 
NSR=12.7% 
Ile(t)ll 2 = 1.664 
NSR=25.7% 
Ile(t)ll2=3.330 
NSR=51.5% 
Ile(t)ll2=4.938 
NSR=76.3% 
Ile(t)ll2=6.660 
NSR= 102.9% 


Table 2.9: AWLS/MFT with randomized X(0)' and Fb = 0.4Hz. 
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Figure 2.3: Normalized Bias and STD plots with randomized ^(0) 
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ness to the randomized initial conditions should not be a surprise if we bear in 
mind the fact that the MFT itself was originated in a way that the estimation 
of initial conditions can be totally avoided through the modulating process. 
Hence, even though almost half of the system response is composed with the 
transient process, the LS/MFT, WLS/MFT and AWLS/MFT algorithms have 
not been thwarted at all. 

2. In this case, PEM also consumed far more computing time than all the MFT 
algorithms combined, partly due to the two more unknowns introduced by the 
initial conditions. 

3. Among the MFT algorithms, the WLS/MFT and AWLS/MFT, again mani- 
fested improvement through a lower bias and standard deviation. 

We have noticed the relatively large bias in the PEM algorithm even in the low noise 
cases from the above simulation studies. This could be attributed to the fact that 
the PEM was developed in a discrete time framework. Therefore, the conversion 
from the discrete time domain to the continuous domain is a must when PEM is 
applied to a continuous system. The transformation used for this converting process 
could contribute to the noticeably larger bias appearing in these simulation studies. 
Another possible cause could be ascribed as the requisite steady state conditions and 
long data ensemble of the PEM was not met in our simulation studies. 
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2.5 Conclusion 


The WLS/MFT and AWLS/MFT algorithms stemming from two different signal 
models have been devised and analyzed in detail. Under different assumptions, 
WLS/MFT is a maximum likelihood estimator and AWLS/MFT is an approximated 
maximum likelihood estimator. When the additive noise in the output is small, then 
the estimated parameters and covariance matrix from AWLS/MFT are fairly close 
to the results of the true maximum likelihood estimate. Lemma 4 has not only fur- 
ther disclosed the insightful relationships among the covariance matrices, but also 
paved the way for the numerical implementation of the WLS/MFT and AWLS/MFT 
algorithms. The recursive banded-sparse matrix inversion scheme in 2.3.3 provided 
a stabler and more efficient method of inverting the covariance matrices. From the 
simulation and comparison studies in section 2.4, the WLS/MFT and AWLS/MFT 
schemes have improved the previous LS/MFT method in both bias and variance, and 
both achieved a smaller variance than the popular PEM algorithm which has the 
worst bias results. Meanwhile, the simulations in 2.4 also show that the initial condi- 
tions have basically no visible affect on the performance of MFT algorithms, which is 
concordant with the theoretical analysis in section 2.1. Again, it has affirmed that the 
MFT method is a potent tool to cope with the identification problems using transient 
I/O data. 
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Chapter 3 


Continuous Time Model 
Reduction Using AWLS/MFT 
Algorithm 


3.1 Overview 

Simplifying a high order or complex model with a lower order model has been deemed 
as one of the most important topics in automatic control, signal processing and other 

engineering and science areas. For many complicated high order models, the reduction 

i 

not only can significantly facilitate their analysis and design, but also makes the digital 
or analog simulation and implementation possible and affordable. In the sense of 
approximation, the lower order models should be able to replicate the time domain, 
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e.g., impulse response or step response, as well as frequency response, e.g., Bode 
diagram, as closely to those of the high order model as possible. From the practical 
point of view, the following situations are very common in demanding simplified 
models: 

1. Given a higher order system, ask reduction. 

2. Given a complicated-looking Bode plot, request a simpler parametrized model. 

3. Given an input/output data pair from an unknown-parameterized higher or 
nonlinear model, demand a model with a specified lower order. 

During the last two decades, many research results on continuous time model re- 
duction have been reported [25] [48] [50] [24] [14] [23]. Most of the earlier work, 
categorically named as classical reduction methods (CRM) [14], has been carried out 
based on classical mathematical approximation theories such as the Pade approxima- 
tion [50] [48], the continued-fraction method, and the time- moment-matching method 
[57]. The essence of CRM schemes is expanding the original system into a Taylor 
series about the origin or the low frequency end, while neglecting the rest of the 
frequency range. This naturally incurs the possibility of low accuracy in the higher 
frequency band and potential loss of stability of the reduced system although the 
original system may have been stable. In order to obtain a stable model, many mod- 
ified reduction schemes, like the stability-criterion and differentiation methods [23] 
can guarantee that the final low order model is stable by allotting some stable poles 
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to the denominator beforehand while letting the numerator be determined by the 
CRM. Unfortunately, in most cases, the accuracy in this kind of algorithm has fallen 
prey to the higher priority of stability. The FF-Pade (Frequency-fitting coupled with 
Pade approximation) method [14] has been proposed to alleviate these drawbacks by 
fitting the mid-frequency range. 

Another class of schemes belonging to the time domain methods [25] computes the 
parameters of a reduced order model so as to minimize a certain criterion function 
characterized by the difference of time domain responses (typically impulse responses) 
to a given driving signal [25]. A reduction algorithm developed by Sakr and Bahgat. 

[47] to obtain an optimal reduced order model for a power plant has been one of the 
examples of this kind of time domain approach. 

Stemming from principal component analysis and singular value decomposition, Balanced- 
Realization [24] has proven to be a seasoned order reduction method in both theory 
and practice; it has been commercialized in the popular MATLAB control toolbox. 

The Balanced- Realization scheme was derived from the “signal injection” viewpoint 
by characterizing the relevant subspaces in terms of responses to injected signals. The 
model reduction is done by eliminating subsystems associated with small singular val- 
ues. 

One common ingredient of all the above schemes is that either high order transfer 
functions or the state space models must be known in advance in order to carry out 
the reduction. This means that they are only used to cope with Situation 1 listed 
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earlier. 


The Modulating Function Technique (MFT) has been no stranger to this field. In [26], 
three nonparametric modulating function frequency matching (MFFM) schemes have 
been initiated and compared with other algorithms. MFFM is a two-step scheme: first 
estimate the frequency response of the system through an nonparametric modulating 
function algorithm [36] from several I/O data pairs; secondly, estimate the parameters 
of a reduced order model by minimizing a frequency matching criterion. Parametric 
LS/MFT has also been resorted to in chemical system reduction [6], though MFFM 
outperformed LS/MFT in [26]. With the more sophisticated AWLS/MFT algorithm 
presented here, better results are expected and will be demonstrated. 


In this Chapter, the AWLS/MFT algorithm stated in the last chapter will be utilized 
to reduce the orders of higher order systems, and the comparison will be made with 
other published results and algorithms, especially with balanced realization, FF-Pade 
and MFFM schemes. As a final example, a 12th order power plant model given by 
Safer and Bahgat [47] is reduced to a 2nd order model and compared with their time- 
domain-based reduction results. The versatility and flexibility of AWLS/MFT will 
be addressed in terms of handling the situations listed above. Some important key 
points of implementation will be mentioned as well. 
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3.2 AWLS /MFT and Order Reduction 


AWLS /MFT has been developed to tackle continuous time parameter identification 
problems based on an I/O data pair. Parameters are obtained through minimizing 
a structure-specified time-domain differential equation model error modulated b} a 
set of known modulating functions. Due to the Fourier type modulating function- 
als we have used, both time and frequency domain information has been naturall} 
concatenated into the joint cost function (2.91). This viewpoint of AWLS/MFT is 
close to what model reduction is trying to accomplish. The essence of the two is so 
indistinguishable that AWLS/MFT should be able to shoulder the mission of model 
reduction. AWLS/M^FT was developed as an I/O-data-pair-oriented scheme, which 
implies that for a given I/O pair of a high order or complex system (the exactly 
parametrized model of this system might not be available), AWLS/MFT can be ap- 
plied to produce a model with a specified lower order. If the transfer function or 
state space representation of a high order system is given and a lower order model is 
desired, the I/O data pair can be acquired by driving the original system with a rich 
persistent signal, typically a Gaussian distributed random sequence, and simulating 
the output with the help of the LSIM routine of MATLAB, so that AWLS/MFT can 
then be used. Many mechanical systems or components are often characterized by 
weird-looking Bode diagrams attached to them, when they are sent out of the man- 
ufacturing or testing sites. If connecting these systems into control loops is desired, 
reduced parametrized models have to be acquired first. In this case, the driving signals 
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could be formed by frequencies coinciding with the modulating function frequencies 
which are decided upon from the Bode diagrams, while accordingly randomizing the 
corresponding phases. One judicious way of choosing the phases would be by using 
so-called low-noise noise [41] which is accomplished by picking phases minimizing the 
fourth moment of the driving sequence. This low-noise noise can smooth out the giant 
peaks which might otherwise drive the system into the nonlinear regions of system 
operation. (This might not be the case in model reduction.) The resulting low-noise 
signal makes the input look like random noise. Gaussian distributed randomized 
phases have been utilized in our simulation studies, and the results have been very 
satisfactory. Another major concern often encountered in model reduction is whether 
or not the reduced order model is stable. This concern can be easily eased with the 
stability-constrained AWLS/MFT algorithm described in the previous chapter which 
automatically locates all the poles of the resulting model in the stable region. 

Technically, the following algorithm-related parameters must be specified before AWLS/MFT 

can be applied 

Fs : Sampling Frequency 
N : I/O Data Length (Number of samples) 
u o : Resolving Frequency = 2n ■ Fs/N 

ub • Modulating Bandwidth (System Bandwidth Covered by Modulating Frequencies) 

M : Maximum Modulating Frequency Index = “integer-part (u>b /<^o)” 

Among the above, u > o is the most intrinsic algorithm-related parameter which de- 
termines the frequency resolution, especially when accuracy is demanded in the low 
frequency range. In some cases in which high accuracy is desired in the middle or 
high frequency range, u> 0 can be set relatively larger so as to alleviate the computation 
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burden while preserving good accuracy. This flexibility is exclusively possessed by 
MFT-based schemes. From the several examples we are going to present here, the 
roles of these parameters will be further illuminated. 


3.3 Comparison with Other Methods 

In order to provide some quantitative measures of goodness in replicating the original 
systems characteristics, the signal-to-error ratios (SER) or S j E are defined, in both 
time and frequency domains respectively, as 

SER< = 20 • logj 0 { , » n time domain t G [0, T s ] (3.1) 

l llM‘j - J 

in frequency domain u> € [uq,^] 

(3.2) 

where || • || 2 denotes the L 2 norm in the appropriate space, and 

h 0 (t) : time response (e.g., step or impulse response) of original system. 

h r (t) : time response (e.g., step or impulse response) of reduced system. 

H 0 (ju>) : frequency response of original system. 

H r (ju) : frequency response of reduced system. 

T s : time interval of interest (roughly the settling time of the system). 

uj\ and u >2 : frequency range of interest. 

Without specific mentioning, all the dB numbers in the graphs of this chapter should 
mean SER or S/E values 1 . 

1 In all six examples, the SER; numbers will be calculated at frequency nodes generated by 
the MATLAB routine LOGSPACE(wi, w 2 , 100), where [u>i,w 2 ] is the same as the graphic range 
appearing in each magnitude or phase plot. Therefore, for the high frequency matching in Example 
1 and 2, those SER numbers are less indicative due to less concern about the low frequency range. 


SER; = 20 • log 


10 


IlgoMlk 

\\H 0 (ju>)-H r (ju)\\ 2 , 
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Example 1: As the first example of our comparison studies, a sixth order low-pass 
system: 

(3 + 2)Ms + 5) 2 -(s+100) 

(s + l) 2 (» + 10) 2 (i + 100) 2 

from [26] is chosen to be reduced to a 3rd order model. Using this system, Co and 

Ydsti [6] compared the FF-Pade algorithm with LS/MFT scheme and found that 

LS/MFT worked better than the FF-Pade method. In Pan’s thesis [26], MFFM gave 

a better fit than the LS/MFT scheme with the following low frequency matching 

model: 




MFFM 


4.5233s 2 -f 30.6739s + 45.2871 
s 3 + 51.2456s 2 + 710.2395s + 453.6771 


(3.4) 


which was obtained through eight-seconds of data at a sampling rate F a = 200Hz, 
i.e., 1600 simulation points. From MATLAB, the reduced 3rd order model using 
Balanced-Realization routine MODRED has been acquired as 


H?- R (s) = 


5.4156s 2 + 30.9466s + 58.5745 
s 3 + 65.4461s 2 + 793.1840s + 585.7454 


(3.5) 


Choosing F s = 64(Hz), u>b = 107r(rad/s) and N = 1024 (resolving frequency: = 

0.1257r(rad/s) 2 ) and using a Gaussian random noise sequence to excite the H\(s), 
2 In this case, the curvature of the low frequency band is slowly-changing, so that wo = 
0.1257r(rad/s) should be fine enough. If uo is further reduced by raising TV, it will not make any 
significant difference. Decreasing u> 0 by lowering the F s is not recommended by and large, due to 
its potential influence on the accuracy in numerical integration. Please refer to Examples 3~6 in 
this Chapter for the cases in which a much finer resolving frequency is required to identify sharply- 
changing peaks and valleys in the Bode diagrams. 
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AWLS/MFT gives the low frequency matching model: 


- AWLS _ 4.7801 5 2 + 28.96235 + 50.4280 ,g ^ 

Hl ~ 7 + 57.8232s 2 + 713.7344s + 502.0671 

The Magnitude, phase and step response plots of the models from Balanced-Realization, 
MFFM and AWLS/MFT algorithms are summarized in Figures 3.1(a), 3.1(b), 3.2 
and 3.3 respectively. In the frequency domain, the three algorithms have very 



10 -2 10" 1 10 ° 10 1 10 2 10 3 10 4 
Frequency (rad/s) 

Figure 3.1(a): Log-scale magnitude plots of the low frequency matching models. 

close SER numbers, though the AWLS/MFT has a slight numerical edge. Graphi- 
cally, especially from the linear-scale magnitude plot and phase plot, the AWLS fits 
best in the low frequency side (about a 2dB lead). From the step response plots, 
Balanced-Realization and AWLS/MFT are better than MFFM while the visual dif- 
ference between Balanced-Realization and AWLS/MFT is negligible. Numerically, 
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Figure 3.1(b): Linear-scale magnitude plots of the low frequency matching models. 



Figure 3.2: Phase plots of the low frequency matching models. 
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Figure 3.3: Step responses of the models, 
the Balanced-Realization has a very small lead (about 0.5dB). 

As mentioned before, MFFT and AWLS are capable of reducing the original system 
to a frequency-range-specified low order model through freely choosing the frequency 
range to match. Pan [26], with this very same example, manifested this kind of 
flexibility and gave a high frequency matching model: 

H MFFM (s) _ 1.4952s 2 + 863.1844s + 5334.4197 

{ s) ~ s 3 + 187.4456s 2 + 11304.2301s + 141250.0129 

With F, = 1024(Hz), u B = 1007r(rad/s), N = 256 (resolving frequency: u> 0 = 

87r(rad/s) 3 ) and a Gaussian random noise sequence exciting the H 1 (s), AWLS/MFT 

3 Like the low frequency matching, the high frequency band also looks “smooth” in this case, and 
the value u>o = 8ir(rad/s) is fine enough. By changing N to make wo relatively coarser or finer, it 
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came up the high frequency (roughly from 87r ~ 1007r(rad/s)) matching model: 

it awls t \ 1.010s 2 + 1003.88025 + 6415.0895 

1 ^ “ s 3 + 214. 2575s 2 + 12816.7076s + 148767.3293 1 j 

In order to see the high frequency matching of the MFFM and AWLS schemes 
clearly, the frequency responses of the two reduced models are drawn together in 
Figures 3.4(a), 3.4(b) and 3.5. Clearly, AWLS/MFT is favored both graph- 



Frequency (rad/s) 


Figure 3.4(a): Log-scale magnitude plots of the high frequency matching models. 


ically and numerically. Another aspect we should notice is that AWLS/MFT used 
only 256 points, equivalent to 250ms of I/O data, and took practically no time to get 
its model, while MFFM used 1600 points which is eight-seconds of data. Therefore, 


will not cause any noticeable change. However, if there is a kink in the high frequency range like 
Example 2 of this Chapter, a finer wq is a requisite. 
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Figure 3.4(b): Linear-scale magnitude plots of the high frequency matching models. 



Figure 3.5: Phase plots of the high frequency matching models. 
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AWLS/MFT can do high frequency fitting more efficiently. 

Example 2: We now consider the following sixth order model [26]: 

4.5s 5 +16.8750s 4 + 1.2474x10V+5.3034x10V+8.0454x10 6 s + 3.6556x10 7 
H ^ s > ~ s 6 + 66.85s 5 + 2778.9s 4 + 7.1963xl0 4 s 3 -|-1.0168xl0 6 s 2 -|-9.8011xl0 6 s-|-4.7306xl0 7 

(3.9) 

Pan again used eight-seconds of data and 1600 data points and by weighting the 
middle frequency range slightly more, came up with the following reduced third order 
model: 

ttMffm i 3.8121s 2 - 59.6845s + 6913.6012 

2 “ s 3 + 42.7534s 2 + 389.3838s + 8388.2031 

This model is basically obtained through a high frequency matching, while he did not 

provide a low frequency match model in this example. From MATLAB, the reduced 

model using Balance-Realization is given by 


(3.10) 


H 2 B ~ R (s) = 


-5.5419s 2 + 170.4683s + 872.5521 
s 3 + 15.4064s 2 -1- 207.4581s + 1129.1082 


(3.11) 


Like the last example, we can still exploit the flexibility of MFT algorithms by fitting 
different frequency ranges. For a low frequency range, setting F s = 256(Hz), u>b = 
27r(rad/s) and N = 1024 (resolving frequency:^ = 0.57r(rad/s)), we have 


fiAWLS 


(s) = 


-4.521 s 2 + 155.179s + 704.284 
s 3 + 14.216s 2 + 196.103s + 914.840 


(3.12) 


If we change the setting to F, = 1024(Hz), ojb — 50?r(rad/s) and N = 1024 (resolving 
frequency:^ = 27r(rad/s)), the high frequency matching model is obtained as 


ijAwls, _ 4.5342 s 2 + 5.1391s + 6581.3344 

2 _ s 3 + 45.0462s 2 + 385.1122s + 6394.8751 


(3.13) 
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Figure 3.6(a): Log-scale magnitude plots of the models. 



Figure 3.6(b): Linear-scale magnitude plots of the models. 
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Figure 3.7: Phase plots of the models. 



Time (s) 


Figure 3.8: Step responses of the models. 
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The performance of these models can be evaluated from Figure 3.6(a), 3.6(b), 3.7 
and 3.8. For the low frequency range in the frequency domain, AWLS/MFT leads 
balanced-realization numerically (18.11/15.91) and the only difference that can be 
observed is that in Figure 3.6(b) Balanced-Realization is off around the peak re- 
gion, while AWLS/MFT coheres to that kink. For the high frequency matching, 
AWLS/MFT does not have a biased peak like MFFM shown in Figure 3.6(b), while 
a magnitude deviation from the original system in the very high frequency range ap- 
peared in Figure 3.6(a). As for that notch concern, AWLS(H)/MFT is much closer 
than the other, as seen in Figure 3.6(a); numerically, AWLS/MFT is better as well. 
In the step response plots, the H? wls (s) has the best numerical result, followed very 
closely by the Balanced-Realization scheme. It is no surprise that both of the high 
frequency matching models, //^ M/ ^ 5 (s) and (s), are not even close to the 

true step response, because after the initial shoot-up the step response is primarily 
determined by low frequency characteristics of the system. Again, AW 1.8/ M FT took 
a much shorter length of data in acquiring the above high frequency matching model. 

Example 3: The following sixth order model was the second example used in [14] 
to compare with other classical Pade approximation methods: 

(1 + 2.0587s)(l + 2.5529s + 5.4342s 2 )(l + 3.2648s + 2.1476s 2 ) 

= (1 + 3.0092s + 0.7970s 2 )(l + 6.8538s + 0.6965s 2 )(l + 0.1394s + 0.6861s 2 ) 

(3.14) 

The FF-Pade scheme reduced the above sixth order system into the following third 
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order model: 


ffFF-Pade^ = 


1 - 1.4257s + 4.3109s 2 
1 + 0.7003s + 0.8613s 2 + 0.0837s 3 


(3.15) 


The MATLAB routine of the Balanced-Realization scheme gave the reduced model 


as: 

B _ R , _ 110.8197s 2 + 26.4163s + 38.2340 
3 " s 3 + 26.7135s 2 + 7.4733s + 38.2340 ^ ‘ ' 

With F, = 256(Hz), u>b = 57r(rad/s) and N = 8192 (resolving frequency: ujq = 

0.06 257r(rad/s)), AWLS/MFT produced the third order model: 


jjAWLS 


(«) 


67.3600s 2 + 6.4604s + 18.8969 
s 3 + 13.7774s 2 + 4.8674s + 18.8969 


(3-17) 


The frequency and time responses are plotted in Figures 3.9(a), 3.9(b), 3.10 and 3.11. 

The AWLS/MFT and Balanced-Realization performed reasonably well, though 
AWLS has shown a slight edge in the frequency domain while Balanced-Realization 
has led in the time domain. The FF-Pade lags far behind numerically and graphically. 
Frankly to say, this is the toughest example for AWLS/MFT we have met in all our 
model reduction studies in the sense of data length required for a very fine resolving 
frequency to identify that narrow valley-peak transition band. 

Example 4: Another middle and high frequency range model used in [14] is the 
sixth order high-pass system: 

1 + 8.8818s + 29.9339s 2 + 67.087s 3 + 80.3787s 4 + 68.6131s 5 
,S ’ ~ 1 + 7.6194s + 21.7611s 2 + 28.4472s 3 + 16.5609s 4 + 3.5338s s + 0.0462s 6 ' 

(3.18) 

Using the middle range frequency fitting and Pade approximation for the lower fre- 
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Figure 3.9(a): Log-scale magnitude plots of the models. 
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Figure 3.9(b): Linear-scale magnitude plots of the models. 
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Figure 3.10: Phase plots of the models. 



Figure 3.11: Step responses of the models. 
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quency band, FF-Pade had the following reduced third order model: 


H FF-Pade^ = 


1 + 2.0098s + 3.7169s 2 
1 + 0.7474s + 0.1898s 2 + 2.4977s 3 


(3.19) 


It should be noted that this reduced model H FF ~ Pade (s) is an unstable system, so 
that it would be unfair to compare its time domain response with the others. The 
Balanced- Realization model reduction scheme laid out its own reduced model as 


H*~ r (s) 


1551.8019s 2 + 554.9813s + 390.9247 
s 3 + 80.3317s 2 + 299.7585s + 390.9247 


(3.20) 


Armed with the F s = 128(Hz), oo B = 27r(rad/s) and N = 8192 (resolving frequency: 
ujo = 0.03127r(rad/s)), AWLS/MFT has reduced the H 4 (s) to 


H. 


AWLS 


(s) = 


1466.9023s 2 + 342.8845s + 355.4733 
s 3 + 75.2876s 2 + 281.7353s + 319.1008 ’ 


(3.21) 


which is a stable system. The frequency domain comparisons are shown in 

Figures 3.12(a), 3.12(b), 3.13. The step responses of the Balanced-Realization and 
AWLS/MFT systems are plotted in Figure 3.14. 

In comparing the frequency domain results for AWLS/MFT and Balanced-Realization 
algorithms, FF-Pade is absolutely in no sense a comparable method. Due to the 
slightly better fitting of AWLS/MFT around the pass-band peak area appearing in 
Figure 3.12(b), AWLS/MFT has a higher SER value. In the time domain response 
of Figure 3.14, both methods give almost a perfect match, while AWLS/MFT has a 
negligibly small numerical lead. 
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Figure 3.13: Phase plots of the models. 



Time(s) 


Figure 3.14: Step responses of the models. 
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Example 5: The last example used in [14] is a complex sixth order system: 


1 + 7.7617s + 13.5756s 2 + 67.6016s 3 + 40.2492s 4 + 144.0994s 5 
1 + 14.8243s + 75.7619s 2 + 163.2959s 3 + 139.3768s 4 + 38.6263s 5 + 3.3282s 6 ' 

(3.22) 


The FF-Pade method gave the reduced fourth order model: 


H[ F - Pade (s) = 


1 + 1.1483s + 4.5589s 2 + 5.0011s 3 
1 + 8.2109s + 5.2508s 2 + 1.4141s 3 + 0.200s 4 ' 


(3.23) 


Again from MATLAB, the Balanced-Realization scheme gave 

b-r, 50.9400s 3 + 68.1574s 2 - 8.3170s + 116.8597 

5 - s 4 + 15.9582s 3 + 56.3846s 1 + 114.1988s + 16.8599' 


(3.24) 


Running AWLS/MFT with F s = 32(Hz), u>b = 27r(rad/s) and N — 4096 (resolving 
frequency: u> 0 = 0.01567r(rad/s)), the following reduced fourth order system has been 
reached: 


H 


AWLS t 


‘iO.OUtfUS TO.L'tlLS -r j.nwo I U.UU4.T : 

= s 4 + 11.5040s 3 + 39.1397s 2 + 39.8916s + 3.1557 




The frequency and time responses of all models are presented in Figures 3.15(a), 3.15(b), 
3.16 and 3.17 respectively. In this example, even though FF-Pade has a relatively 
good fit in the notch part of Figure 3.15(a), the overall AWLS/MFT and Balanced- 
Realization are still better as shown in Figure 3.15(b) and 3.16, while AWLS/MFT 
has the closest frequency fitting, especially in the peak area, numerically and graph- 
ically. In the time domain, except for the laxge over-shoot of FF-Pade at an early 
stage, they all agree well with the true step response, though the Balanced-Realization 
has a slight lead over AWLS/MFT numerically. 
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Figure 3.15(a): Log-scale magnitude plots of the models. 
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Figure 3.15(b): Linear-scale magnitude plots of the models 






Figure 3.17: Step responses of the models. 
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Example 6: As the final example of our model reduction studies, we shall reduce 
a twelfth order power plant system arisen in a power generating station [47] to a 
non-strictly-proper second order model and compare it with the results stated in [47]. 
Sakr and Bahgat gave the following transfer function bridging between the output 

power and the input disturbance: 

684s 11 - 93484.9s 10 - 1.76 x 10 6 s 9 - 3.25 x 10V - 3.336 x 10V 
He ^ = s 12 + 16.9s 11 + 335.88s 10 + 3694.9s 9 + 30709s 8 + 199002s 7 + 928769s 6 


-1.944 x 10 9 s 6 - 7.079 x 10 9 s 5 - 1.697 x 10 lo s 4 - 2.71 x 10 lo s 3 
+3.02 x 10 6 s 5 + 6.81 x 10 6 s 4 + 1.047 x 10 7 s 3 + 1.056 x 10 7 s 2 


-2.80 x 10 lo s 2 - 1.703 x 10 lo s - 4.926 x 10 9 og . 

+6.344 x 10 6 s + 1.7121 x 10 6 

Based on minimizing the cost function constructed with the time domain responses, 
the parameters of a non-strictly proper 2nd order system is optimally computed and 
this method will be denoted as the S-B scheme in our following discussion. In [47], 
the S-B algorithm with steady state constraint seems to give much better results than 
the one without the steady state constraint, so that only the reduced model with the 
steady constraint will be compared with the AWLS/MFT here. With the steady state 
constraint, S-B has reduced He(s) to 


Ht B (s) = 


6365.3s - 142182.23 
s 2 + 1.475s + 52.311 + 


(3.27) 


Using BALREAL and MODRED routines of MATLAB, Balanced-Realization gives 


Ht R (s) = 


-266.884s 2 + 1112.084s - 142495.032 
s 2 + 0.5246s + 49.5261 


(3.28) 
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Setting N = 2048, u>b = 2x and Fs = 50Hz (resolving frequency: u>o — 0.0487r(rad/s)), 
AWLS/MFT has the reduced model: 


68.477s 2 + 532.948s - 125098.5 
s 2 + 0.5235s + 49.3780 


The frequency and time responses of the reduced models are plotted in Figures 3.18(a), 
3.18(b), 3.19 and 3.20. 
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Figure 3.18(a): Log-scale magnitude plots of the models. 


In the time domain, Figure 3.20, where S-B has been the derived, the AWLS/MFT 
and Balanced- Realization have accurately located the impulse response of the original 
system, while the impulse response from the S-B model is dying out too fast and its 
discrepancy from the true system is obvious, both visually and numerically. From 
the frequency response plots, Figures 3.18(a), 3.18(b) and 3.19, the AWLS/MFT and 
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Figure 3.18(b): Linear-scale magnitude plots of the models. 
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Figure 3.20: Impulse responses of the models. 

Balanced-Realization also score significantly higher than the S-B scheme. Between 
AWLS/MFT and Balanced-Realization, the AWLS/MFT has better SER numbers 
in both the frequency and time domains, though the graphical difference is almost 
invisible. Overall, the poor performance of the S-B method is mirrored in front of the 
AWLS/MFT and Balanced-Realization schemes. 


3.4 Concluding Remarks 

From the first two examples, AWLS/MFT did show the promising improvement over 
the MFFM schemes in both accuracy and data source efficiency, especially when high 
frequency matching is desired. The FF-Pade method has been compared with the 
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AWLS/MFT and Balanced- Realization algorithms in three examples and has failed 
to prove that it is a competitive alternative to these two algorithms. Meanwhile, the 
FF-Pade still faces the possibility of the resulting model being unstable. In the final 
example, the S-B method was almost a fiasco relative to the results of the AWLS/MFT 
and Balanced- Realization algorithms. The Balanced-Realization has been utilized in 
all the examples and it has held quite a good stance in both the time and frequency 
domains. It is safe to say that the AWLS/MFT algorithm is working at least as well 
as the Balanced- Realization scheme due to the fact that in the frequency domain 
the AWLS/MFT has more or less an edge (larger SER numbers) over the Balanced- 
Realization, while in the time domain they both share a tiny numerical lead one 
way or another. One interesting point that should be noted regarding the Balanced- 
Realization scheme is that model reduction is just a very narrow application of the 
Balanced- Realization technique. Quoting from [24]: 

“the relationship between general model reduction and reduction by sub- 
system elimination is not well understood”. 

This somehow is reminiscent of the vagueness of applying AWLS/MFT to model 
reduction problems. Even though it was developed as a parameter identification al- 
gorithm, the AWLS /MFT has shown in the above examples that it can indeed fulfill 
the mission of model reduction. The Balanced-Realization method takes almost no 
time to perform the model reduction, while the overall time (including synthesizing) 
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of AWLS/MFT is ranging from no time to 6 seconds in above examples 4 , which does 
not pose any formidable threat to its practical usage. The limitations of the Balanced- 
Realization reduction scheme compared with the MFT methods include: (i) although 
its results fit the low frequency range well, it lacks the flexibility of choosing the fre- 
quency range to match models as is possible in the AWLS and MFFM methods; this 
is demonstrated in the first two examples, and (ii) one premise of applying Balanced- 
Realization is that the original exact state space or transfer function model must be 
known beforehand. In this sense, AWLS/MFT algorithm has injected some sort of 
versatility into the model reduction technique. We have successfully tried to use the 
Bode diagrams of the above example models instead of analytical high order trans- 
fer functions as the start point to approximate these Bode diagrams with specified 
(low) order models and found that the results are almost identical with those started 
with high order transfer functions. This underscores the flexibility and versatility 
of AWLS/MFT in model reduction as one of the major features that could make it 
stand out among its peers. 


4 This is directly related to the resolving frequency used in reduction. Usually the high frequency 
matching takes much less time than the low frequency matching reduction. 
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Chapter 4 


Parameter Estimation of MIMO 
Systems Without Pole Constraints 

4.1 Introduction 

Consider a multi-input and multi-output (MIMO) continuous time system given by 
the following transfer function form: 

y(s) = H(s)u(s) (4.1) 

where u(s) is a pi x 1 transformed input vector, y(s) is a p 2 x 1 transformed output 
vector, H(s ) is a pi x p 2 transfer function matrix with ( k,q) th element 

**■« = (4-2) 

and ( Bk q {s), y4fc g (s)) are coprime polynomials in s. Similar to the notation of the SISO 
cases in Chapter 2, denote 0k q as the parameter vector formed by the coefficients of 
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the i4i,(s) and Bk q {s) polynomials, e.g., for an n th order polynomial: 


— S n + Gfcj.js" 1 + • • • + Ofcg.n 


and a (n — l) th degree polynomial: 


Bkq(s) — bkq t lS n 1 + - ‘ ' + b kq,n 


then 


&kq = 


( a kq,l \ 

~ a kq, n 

b kq,l 

V b kq,n ) 


(4.3) 


Our goal in this Chapter is a systematic procedure of using or extending the modulat- 
ing function technique (MFT) to estimate the totality of parameter vectors Okq based 
on the input-output data pair {u(t),y(f)}, t e [0,T]. As discussed previously, the 
order n is presumed chosen beforehand. In practice, this order may be determined by 
re-solving the identification for increasing orders, starting from an initial value, until 
the residuals and/or SER’s values are sufficiently small and/or large. 
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4.2 Parameter Identification of MIMO Systems 
Using AWLS/MFT 

4.2.1 Regression Form of Modulated MIMO Systems 

Considering the k th output ^(s), we have 

" w “58i)“* w - (4 - 4) 

Denote A*(s) as the least common multiple of {/4*i(.s), Ak Pi (s)} and 

(4 - 5) 

Then a modulatable higher order differential operator form of equation (4.4) can be 
written as 

Mp)yk(t) = Y B kg (p)u q (t) (4.6) 

?=1 

where argument p represents the differential operator p = Using a similar notation 
and adding the error terms to fit in the identification framework, system (4.1) is then 
represented in the modulatable differential form: 

f A.lpJy.W = E ", Bi,(p) U ,(!) + e,(f) 

I . : _ (4 ' 7) 

V A n (p)y P2 {t) = Ylq=l Bp 2 q{p)u q {t) e P2 {t) 

Letting {n* ; k = 1,2,- - • , P 2 } be the corresponding orders of the polynomial set 

{j4fc(p) ; k = 1,2, • • • ,$> 2 } and {@k ; k = 1,2, • • • ,p 2 } the parameter vectors compris- 
ing the coefficients of the polynomials {Ak(s), Bk q (s) •, q = 1,2, as in (4.3), 

k = 1, • ■ • , p 2 i apply the complex form of the modulating function set (2.6) of orders 
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{n* ; k = 1,2, • • • ,p 2 } successively to the above equations. Notice that these models 
axe overparametrized (in general) due to the cross products of the underlying polyno- 
mials comprising the h kq (s). Similar to the SISO case in Chapter 2 (see (2.12), (2.13) 
and (2.14)), (4.7) could be modulated into the following regression set: 


Y x 


fi^i + Cl 


(4.8) 


l ^P2 “ ^P2^P2 “I" ^>2 

If further we denote by {M k ; k= 1, 2, • • • ,p 2 } the highest modulating frequency index, 
then each error vector l k in (4.8) is of dimension (M k + 1) x 1. 

Remark: Here we need to indicate that this problem will be solved as p 2 distinct 
2-stage problems: (i) obtain the {0 k ; k=l,2,- • < ,p 2 ) for the over-parametrized mod- 
els, (ii) reduce each over-parameterized model to the original 0 kq through the model 
reduction scheme discussed in Chapter 3. 


4.2.2 Joint Likelihood Cost Function 


For simplicity, we first look at the cases when {efc,f k = 1,2, • • • ,p 2 } are real 1 . 
If the error sequences {e k ~ ^(0, W k ) ; k = 1,2, • ,p 2 } are mutually independent of 

each other, then the joint likelihood density function can be written as 


1 


. exp {-\ilW k 'h} 


P2 _ P2 

■ J (2 ~ ~ w*» 


^his results by using only the real or imaginary part of the modulating function set to modulate 
( 4 . 7 ) 
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1 


- - exp{ — — J{0\, • • • , 0 P2 )}-(4.9) 

(2 7r )(M 1 +-+M P2 + P2 )/2| W/i |i / 2...| l y i)2 | 1 /2 2 


The log-likelihood C(8 U - • • ,0 P2 ) can be written as 


V 2 


C(8 u -‘-,e p2 ) = ln{ n p(e fc |6> fc )} 

k = 1 


1 


= J(«i, 


(4.10) 


where J(#i, • • • ,0 P2 ) is a quadratic function 


P2 


j(0x, ■■■,«„) = E (n - r^fwrin - w*) 


*=1 

> 0 


(4.11) 


and 

1 P2 1 P2 

K(Wi, • * * , W P2 ) = -i • (p 2 + £ Mfc) • ln(27r) - - £ In |W*|. (4.12) 

1 k = i z fc=i 

If the {W* ; k = 1, 2, • • • ,p 2 } are known, then maximizing the log-likelihood function 
(4.10) is equivalent to minimizing the quadratic function J{9\ , , 0 P2 ), and in this 

case the minimizing parameter set: 

(#i, • • • ,# P2 ) = arg ,min J(#i, • ■ • , 0 P2 ) (4.13) 

will lead to the maximum likelihood estimate of the parameter set {0* ; k = 1 , 2, • • • , P 2 } • 
Regarding the two signal models, i.e., the equation error model and the measurement 
noise signal model, a few comments could be made here: 

1. For the measurement noise signal model without input additive noise, i.e., 
{v q (t) = 0; q = 1,2, • • • ,pi} in Figure 4.1, the condition that the {e* ~ -A/"(0, IF*); 
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Figure 4.1: MIMO measurement noise signal model. 

k — 1,2, •• • ,p 2 } 2 are mutually independent is automatically met; hence, (4.10) 
is the desired joint cost function. Using Lemma 4 of Section 2.3 which guaran- 
tees that the real and imaginary parts of e* are uncorrelated, the terms in (4.11) 
can be directly adjusted to the combined forms (see Equations (2.95)~(2.97)) as 
discussed in the implementation Section 2.3. Neglecting the function k(Wi, • • • , W P2 ) 
for computational expediency, the result from (4.13) is, again, just an approxi- 
mated maximum likelihood estimate. In reality, the additive noise level in the 
inputs of many physical systems is significantly lower than that in the outputs. 

One typical example is the flight data of an F-18 jet aircraft where the inputs 
include longitudinal pilot stick deflection, lateral pilot stick deflection, horizon- 

similar to Equation (2.45) in the SISO measurement noise signal model, all terms involving the 
{«,(<); g=l,2, -.pi} vanish. 
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tal rudder deflection, horizontal tail deflection and left/right aileron deflections; 
these are basically free from noise, but the outputs are quite contaminated. 
Therefore, the focus of our studies will be mainly devoted to this case. 

2. When the input additive noises in the measurement noise signal model of 
Figure 4.1 are non-zero, i.e., ^ 0 ,q = 1,2, • • • ,/>i}, the error sequences 

{e k ; k = 1 , 2 , • • • , P 2 } are generally correlated and, in this case, it is no longer 
appropriate to claim that the result is still an approximated maximum likelihood 
estimate. Nevertheless, the estimate from (4.13) might still be acceptable. 

Heretofore, (4.11) will be utilized as a general joint cost function for the MIMO 
parameter identification problem. As discussed in relation to Lemma 2 and 3 of 
Chapter 2, the covariance matrices usually can be decomposed as 

W k = a\W_k k= 1,2, •• • ,p2 (4-14) 

where {a k ; k = 1, 2, • ■ • , p 2 } are the variances of the equation error noise or additive 
output measurement noise, and { VF* ; k = 1,2,--- ,^>2} are explicitly only related 
to the binomial coefficients, unknown parameters and modulating frequencies (see 
Section 2.2.2 and 2.2.3). We assume that the forms or shapes of the probability density 
functions of the noises are given. In most practical problems, neither the variance of 
the additive noises nor their probability density functions are known exactly. Usually 
an assumption can be made about the distribution of noises, but not their variances. 
In this case, introduce non-negative constants {v k ; k = 1, 2, • ♦ • , P 2 } and modify (4.11) 
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to 


J{h, • • • , <y = £ v*(n - f (4.i5) 

k = i 

where the unknown {o* ; A: = 1 , 2, - - - ,p 2 } are absorbed into the {i/* ; k= 1, 2, • • • ,p 2 }. 
These will be determined along with the parameters in the algorithm to be proposed 
below. 


4.2.3 Decomposability of MIMO Into MISO Models 

As indicated in the title of this Chapter, our discussion will be carried out only under 
the assumption that {Ok ; k = 1,2, •••,p 2 } are independent one to another and no 
constraints are attached to them. From the necessary condition of minimizing (4.15): 

^W k — : - f A) = 0 (416) 

we obtain 

h = (Tlw:'r t )-'t T t w;'Y t . (4.17) 

This result shows that the fcth estimated parameter vector Ok is only determined 
by the output of the kth channel of the MIMO system and by the totality of the 
inputs {u,(f) ; 9 = 1,2, ••• ,pi}. The important conclusion that can be drawn from the 
above is that the parameter estimation problem for a pi-input-and-p 2 -output (MIMO) 
system can be decomposed into a total of p 2 independent MISO sub-problems, if the 
joint cost function (4.15) is employed. 

Under the assumption that the output additive noises {«*(<) ; k= 1,2, • • • ,p 2 } are mu- 
tually independent white Gaussian processes and the input additive noises {u 9 (t) ; q — 
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1,2, * * * , Pi } are of zero variance (or significantly smaller than the variances of the 
output noises), the AWLS algorithm derived in Chapter 2 can be almost directly 
applied to obtain an approximate maximum likelihood estimate Ok for the fcth MISO 
sub-system {Ak(s), Bk q (s) ] q= 1,2, ■•*,pi}. Moreover, only those elements in the 
parameter vector Ok related to the denominator polynomial Afc(s) will be involved in 
the adaptive iterative procedure. 

4.2.4 Reduce from {A*(s), Bk q {s)} to {A* g (s), Bk q (s)}; q = 1, • • • ,Pi 

As mentioned earlier, our ultimate goal is to estimate the coefficients {0k q ; q = 

1,2, • • • ,pi} of the original coprime pairs {Afc,(s), Bk q (s) ; 9 = 1,2, • • ■ ,pi} which are 
not directly modulatable; instead, we have formulated a joint cost function for es- 
timating the Ok for higher order modulatable pairs {i4jt(s), Bk q (s ) ; q = 1,2, • • • ,Pi}. 
Therefore, a model reduction problem will be confronted in the process of going 
from ; q = 1,2, • • • ,pi} to the lower order pairs {A kq {s), B kq {s) ; q = 

1,2, • • • ,pi}. Even though there is a little bit of reluctance to admit it, this actually 
is the very reason that triggered our studies on model reduction problems summa- 
rized in Chapter 3. From the six examples in Chapter 3, AWLS/MFT is known 
to be successful in reducing some weird-looking high order systems into lower order 
models with accuracy at least equal to the seasoned Balanced- Realization model re- 
duction scheme. As another typical application, the model reduction required from 
{A k (s),Bk q (s) ; 9 = 1, 2, • * • , pi } back to {A kq (s), B kq {s) ; q= 1,2, • • • ,pi} will be con- 
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ducted using the AWLS/MFT schemes in our upcoming numerical simulation studies. 


4.2.5 Overall Procedure of Estimating MIMO Systems 

After these discussions, we propose the following procedure of using the modulating 
function technique (MFT) to estimate the parameters of a MIMO system, provided 
the system structure of each h kq (s) be known: 

1. In regards to the unmodulatable original MIMO system model (4.4) with its 
polynomials 

{A kq (s),B kq {s ) ; 9 = 1,2, k- 1,2, • • • ,&} 

focus on the higher order modulatable MIMO system model ( 4.6) with its 
polynomials 

{A k (s),B kq (s)-, 9=1,2, • • • ,pi’, *=1,2, ■ • ■ ,p 2 } 

which are obtained by multiplying both sides of (4.4) with the least common 
multiple polynomial {A*(s) ; k= 1, 2, • • • ,p 2 }- 

2. Select a set of nonnegative weights {i/ k ; k = 1,2, •• • , P2} , and utilize the joint 
cost function (4.15) thereby decomposing the MIMO system 

{A*(s),i?* 9 (s); 9 = 1,2, ••• ,pi ; *=1,2, • • • ,p 2 } 

into p 2 sets of MISO sub-systems {A*(s), B kq (s ) ; 9= 1, 2, • • • ,pi}. 

3. Initialize the index k = 1. 
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4. Apply AWLS/MFT to the Arth higher order sub-system {A*(s), Bk q (s ) ; q = 
1,2, • • • ,pi} and get its approximate maximum likelihood estimate Ok- 

5. Using the AWLS/MFT model reduction technique of Chapter 3, convert Ok back 
into {Okg ; <jr = l,2,* •• ,pi} with following steps: 

(a) Inject the q th input u q (t) into the estimated q th higher order model h kq (s) 
with the MATLAB simulation routine to obtain yk q {t), be., 

y kq = LSIM(B fc „ ,«,(*), T s ,), 

Hence, we have acquired an I/O data pair {u q {t),yk q (t)} for the SISO 
system hk q (s). 

(b) Estimate 0k q for the pair (Ajt 9 (s), Bk q {s)) using AWLS/MFT and the I/O 
data pair {u q {t),yk q {t)} from (a). Keep the original algorithm related 
parameters (woi^b) as used in Step 4. 

6. k = k + 1; if k > p 2 , stop; otherwise go back to step 4. 


4.3 Numerical Simulation Results 

Due to the decomposability under the joint cost function J(0\, • • • , 0 P2 ), it is only nec- 
essary to undertake numerical experiments with MISO systems. Meanwhile for con- 
venience, the index k will be dropped from our previous notations for the MISO sub- 
systems, i.e., modulatable {A(s), B q (s ) ; q= 1, 2, • • • ,pi} or unmodulatable {h q (s ) ; q- 
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1,2, • • • , pi } will be sufficient to represent a MISO system. But first, some potential 
problems we might face will be stated. 

4.3.1 Several Possible Combinations of {A(s), B q (s)} 

Questions will naturally arise like whether or not large discrepancies in gain and 
frequency bandwidths among the {h q (s); q = 1,2, • • • ,pi} will cause some kind of 
frenzy in terms of accuracy. For the modulating function technique, it could be a 
major concern that choosing the maximum modulating frequency index M based on 
the maximum bandwidth among {h q (s) ; q = 1,2, • • • ,pi} might degrade the estimate 
for those A,(s)’s with relatively narrow bandwidths. When the only available output 
is corrupted with additive noise, the impact on the &,($)’ s with smaller gains could be 
devastating. In order to explore these potential difficult combinations, five 2-input- 
and-single-output systems { h q (s ) ; <7=1,2} will be simulated in our numerical studies. 
These are configured as follows: 

1. Roughly the same gains and frequency bandwidths between hi(s ) and /i 2 (s); 

2. Roughly the same gains but different frequency bandwidths between hi(s) and 

h 2 (s)] 

3. Roughly the same frequency bandwidths but different gains between (s) and 

M 5 ); 

4. hi(s) has a higher frequency bandwidth and a higher gain than /i 2 (s). 
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5. hi(s) has a higher frequency bandwidth but a lower gain than h 2 (s). 

4.3.2 Setup of Numerical Simulations 

Two hundred Monte Carlo runs are going to be carried out at each of several noise 
levels for each of the five systems mentioned above. The driving inputs ui(t) and 
u 2 (t) are two Gaussian random sequences generated through the MATLAB routine 
RANDNQ which is basically a white noise generator. Using the linear system simu- 
lation routine LSIM(), the outputs of hi(s) and h 2 (s) can be simulated symbolically 
as 


y i(0 

= LSIM (fc 1 (a),ii 1 (0,<); 

(4.18) 

V2{t) 

= LSM(h 2 (s),u 2 (t),ty, 

(4.19) 


Then a single output y(t) available for identification purposes is synthesized through 
superposing yi{t) and y 2 (t) as 

y{t ) = y\(t) + V 2 (t) +n{t) 

' V / 

y(0 

= y(t) + n(t) (4.20) 

where y(t) represents the noise-free (ideal) single output and n(t ) is the additive noise 
formed as 

n(f) = K • STD(yi(i)) • RANDN(); (4.21) 

where the MATLAB routine STD() returns the standard deviation of y\(t) and the 
constant /C controls the amplitude of n(t). In order to see the affect of additive noise 
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on each subsystem, two noise-to-signal ratios (NSR) are defined as 


NSR, = • 100%, * = 1,2. (4. 

When K in (4.21) is zero, which means no additive noise to the single output y(t), 
the two drivings Ui(t) and u 2 (t) will be realized independently for each Monte Carlo 
run. This kind of simulation will evaluate the performance of the algorithm in rather 
ideal situations. But another approach will be adopted here if K- ^ 0. Thus, n(t ) will 
be realized independently in each Monte Carlo run, while the two inputs Ui(t) and 
u 2 (t), generated independently by RANDNQ beforehand, will remain intact during 
the whole ensemble of Monte Carlo simulations. In this way, the noise levels NSRi 
and NSR 2 on yi(t) and j/ 2 (0 will remain basically unchanged during the 200 Monte 
Carlo runs. This arrangement not only facilitates our focus on the “reactions” of 
hi(s ) and ft 2 (s), but also makes some cross comparisons of the five examples possible. 
Like before, the time domain performance of the estimated models can be evaluated 
by the signal-to-error ratio (SER) which is specified as 



SER = 20 • log 


10 


f l|y(*)lla ] 

llltfO-iWIU 


(4.23) 


where y(t ) is the simulated output based on the estimated parameters driven by 
iij(f) and u 2 (t). Note that each SER number will be a random number due to the 
randomness of n(t ) so that its mean and standard deviation will be given as a pair 
of numbers (ttiser, &ser) in our final numerical tables as determined from the 200 
Monte Carlo runs. Let Qk be the kth element of the estimated parameter vector 6 and 
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Q° k as the kth element of the true parameter vector 6°. Two additional normalized 
quantities, 

normalized bias NB(%): NB = |— | • 100% 

9 k 

and 

normalized standard deviation NSTD(%): NSTD = |— • 100% 

9k 

which characterize the accuracy in parameter space, will be plotted and compared, 
besides tabulating the mean (gk) and variance (<7^) for each estimated parameter Qk- 
For all five systems, the running parameters of MFT are chosen as follows: 

• sampling frequency: f s = 25.6Hz 

• number of data points: N = 512; time interval: T = 20s 

• resolving frequency: f 0 = 0.05Hz or ui 0 = 0.314(rad/s) 

• modulating frequency bandwidth: fs = 0.5Hz or ub = 3.142(rad/s) 

• maximum modulating frequency index: M — 10 


4.3.3 Numerical Experiments 

Example 1: Consider the two input and single output system 


y( s ) 


hi(s)ui(s) + h 2 (s)u 2 {s) 

s 2 + 3s + 8 Ul ^ + s 2 + 4s + 10 


U 2 (s). 


(4.24) 
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The Bode diagrams of /ii(s) and are plotted in Figure 4.2 from which we can 

see that hi(s) and have rather similar bandwidths and gains. Two hundred 

Monte Carlo run results are listed in Table 4.1 from which we can see that if there 
is no additive noise, the algorithm can give an almost perfect estimate, and in this 
example NSRi and NSR 2 are roughly in the same range in each of the different 
noise levels. The normalized bias (NB) and normalized standard deviation (NSTD) 



Figure 4.2: Bode diagram of /i x (s ) and in Example 1 


computed from the table are plotted in Figures 4.3 and 4.4 respectively from which 
it seems that the influence of the additive noises on hi(s) and /i 2 (<s) are relatively the 
same, though hi(s) has a slightly smaller bias. This case will be used as a reference 
in the following comparisons. 
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NSTD(%) NSTD(%) nB(%) NB(%) 


parameter: 5 




parameter; 3 




parameter. 8 




J 


1 

N 


Figure 4.3: Normalized bias vs. additive noise level NSRj 


parameter: 5 




parameter: 3 




parameter 8 




j 


T 
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Figure 4.4: Normalized standard deviation vs. additive noise level NSRi 
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hj(s) 

h 2 (s) 

■ 

m SER 

true 0 

5 

3 

8 

6 

4 

10 



m 

4.9957 

3.0022 

7.9975 

5.9965 

4.0047 

9.9986 

0% 

57.92dB 

CX 2 

0.00001 

0.00000 

0.00001 

0.00002 

0.00001 

0.00002 

0% 

1.567dB 

m 

5.0033 

3.0078 

7.9983 

6.0045 

4.0179 

9.9966 

5.22% 

42.57dB 

a 2 

0.00238 

0.00122 

0.00341 

0.00744 

0.00385 

0.01615 

4.37% 

2.912dB 

m 

5.0581 

3.0245 

8.0606 

5.8998 

3.9530 

9.8033 

9.99% 

35.17dB 

a 2 

0.01723 

0.00770 

0.02462 

0.05194 

0.02939 

0.13041 

8.367% 

3.630dB 

m 

5.1770 

3.0747 

8.2181 

5.5613 

3.6837 

9.3342 

19.48% 

28.1 ldB 

a 2 

0.08889 

0.04878 

0.12089 

0.30429 

0.15721 

036263 


4.345dB 


m: mean value, a 2 : variance, NSR: noise to signal ratio, SER: signal to error ratio. 


Table 4.1: Numerical results of 200 Monte Carlo runs for Example 1 
Example 2: As the second example, consider the two input and single output system 


y( s ) 


^l( 5 ) u l( 5 ) 4" ^2{ s ) u 2{ s ) 

5 , . . 3 


s 2 + 3s + 8 


Ui(s) + 


s 2 + 4s + 10 


u 2 (s). 


(4.25) 


The Bode diagrams of &i(s) and /i 2 (s) are plotted in Figure 4.5 from which we see 
that the gain of h 2 (s) has been reduced from 0.6 to 0.3 while their bandwidths are 
still roughly the same. The results of 200 Monte Carlo runs are listed in Table 4.2, 
and the normalized bias (NB) and normalized standard deviation (NSTD) computed 
from the table are plotted in Figures 4.6 and 4.7 respectively. Bearing in mind the 
corresponding relationship between NSRi and NSR 2 at each additive noise level in 
this case 


relationship between NSRi and NSR 2 

noise 

level 1 

level 2 

level 3 

level 4 

NSRi 

0% 

5.03% 

9.62% 

21.31% 

nsr 2 

0% 

8.04% 

16.10% 

35.68% 
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m: mean value, a 2 : variance, NSR: noise to signal ratio, SER: signal to error ratio. 


Table 4.2: Numerical results of 200 Monte Carlo runs for Example 2 
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NSTD(%) NSTD(%) NB(%) NB(%) 


parameter: 5 



parameter: 3 



parameter: 3 




parameter: 8 




Figure 4.6: Normalized bias vs. additive noise level NSRi 


parameter: 5 




parameter: 3 




NSR1<%) 


parameter: 8 




Figure 4.7: Normalized standard deviation vs. additive noise level NSRi 
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and comparing with Example 1, we have the following observations: (i) Without 
additive noise in y(t ), both /ti(s) and h 2 (s ) can be accurately identified, (ii) The 
NSTD numbers and NB numbers for h 2 (s) have increased by 1% ~ 2% and 3% ~ 4% 
respectively at each noise level, which could be ascribed to the relative larger NSR 2 , 
(iii) Note that the corresponding NSTD and NB numbers for h x (s) have decreased 
by roughly 1% ~ 2% at each noise level. These inverse trends for Ms) and M 5 ) 
exaggerate the contrast in the NB and NSTD plots relative to those in Example 1. 
Example 3: As the third example, we consider the system 


y( s ) 


h\{s)ui(s) + h 2 {s)u 2 (s) 


s 2 + 3s + 8 


Ms) + 


s 2 + 2s + 1 


Ms). 


(4.26) 


The Bode diagrams of h x {s) and Ms) are plotted in Figure 4.8 from which we see 



h'i(s) | 

h 2 (s) 

““““ 

KISH 

ttlsER 

true 0 

5 

3 

8 

0.6 

2 

1 

| 

&SER 

m 

4.9940 

3.0014 

7.9967 

0.5998 

1.9997 

0.99% 

0% 

57.74dB 

CP 

0.00001 

0.00001 

0.00001 

0.00000 

0.00002 

0.00000 

0 % 

2.07 OdB 

m 

4.9970 

3.0051 

7.9911 

0.5976 

1.9928 

0.9%7 

5.26% 

43.34dB 

CP 1 

0.00120 

0.00070 

0.00172 

0.00017 

0.00190 

0.00036 

9.04% 

2.754dB 

m 

5.0039 

3.0070 

7.9761 

0.5893 

1.9654 

0.9859 

10.02% 

37.25dB 

CP 

0.00479 

0.00282^ 

0.00703 

0.00072 

0.00787 

0.00144 

17.21% 

2.823dB 

m 

4.9973 

2.9924 

7.9018 

05575 

1.8621 

0.9468 

19.11% 

30.7 6dB 

CP 

0.01502 

0.00918 

0.03589 1 

0.00293 

0.03144 

0.00527 

32.83% 

2.97 8dB 


m: mean value, a 2 : variance, NSR: noise to signal ratio. SER: signal to error ratio. 


Table 4.3: Numerical results of 200 Monte Carlo runs for Example 3 


that Ms) has a noticeably lower frequency bandwidth while they have roughly the 
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Figure 4.8: Bode diagram of /ii(s) and fi 2 (s) in Example 3 
same steady state gains. The results of 200 Monte Carlo runs are listed in Table 4.3, 
and the normalized bias (NB) and standard deviations (NSTD) computed from the 
table are plotted in Figures 4.9 and 4.10 respectively. Again, an almost flawless 
estimate has been obtained when NSR, = 0. Because of the narrower bandwidth, the 
mapping relationship between NSRi and NSR 2 is changed to the following: 


relationship between NSRi and NSR 2 

noise 



level 3 

level 4 

NSRi 



10.02% 

19.11% 

nsr 2 

0% 

9.04% 

17.21% 

32.83% 


Like Example 2, the noise level imposed on A 2 (s) is almost doubled. But unlike 
Example 2, the following two astonishing observations with respect to Example 1 can 
be made from Figures 4.9 and 4.10: (i) Contrary to intuition, NSTD and NB numbers 
of h, 2 (s) are very similar to and even slightly smaller than those in Example 1. This 
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parameter: 0.6 



parameter: 3 
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Figure 4.9: Normalized bicLS vs. additive noise level NSRi 


parameter: 5 




parameter: 3 




parameter: 8 
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j 

10 


Figure 4.10: Normalized standard deviation vs. additive noise level NSRi 
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implies that narrowing the bandwidth of system h 2 (s ) does not worsen the quality 
of the estimate for h 2 (s) at all, though this narrowing is equivalent to increasing 
the additive noise to it; (ii) For system ^i(s), its NB values has been decreased 
dramatically and its NSTD numbers are also down roughly to half of that in Example 
1. Another important phenomenon in our Monte Carlo simulation processes is that 
the iteration steps required for convergence of the AWLS/MFT algorithm in this case 
are noticeably fewer than those cases in Example 1 and 2 where both bandwidths are 
close to each other. Apparently the AWLS/MFT is responsible for all of these. 
Example 4: As the fourth example, we consider the system 


y ( s ) 


Ms)ui(s) -I- h 2 (s)u 2 {s) 

5 , , .3 


s 2 + 3s + 8 


Ul(s) + 


s 2 + 2s + 1 


u 2 (s). 


(4.27) 


The Bode diagrams of hi(s) and h 2 (s) are plotted in Figure 4.11 


from which we can 


r~ 

hi(s) 


h 2 (s) 

BE M 

rrisER 

true 6 

5 

3 

8 

0.3 

2 

1 

m 

C*SER 

m 

4.9943 

3.0018 

7.9970 

0.2999 

1.9998 

1.0000 

0% 

57.1 ldB 

(J 2 

0.00001 

0.00000 

0.00001 

0.00000 

0.00006 

0.00001 

0% 

1.986dB 

m 

4.9984 

3.0047 

7.9895 

0.2955 

1.9695 

0.9889 

5.03% 

43.5 2dB 

CP 

0.00078 

0.00049 

0.00158 

0.00014 

0.00649 

0.00110 

17.30% 

2.729dB 

m 

5.0072 

3.0065 

7.9655 

02822 

1.8799 

0.9558 

10.10% 

36.65dB 

O 2 

0.00393 

0.00268 

0.00693 

0.00081 

0.03339 

0.00586 

34.72% 

3.322dB 

m 

4.9888 

2.9920 

7.8795 

0.2504 

1.6642 

0.8847 

20.04% 

29.41dB 

G 2 

0.02164 

0.01136 

0.04158 

0.00397 

0.16048 

0.02595 

68.85% 

3.938dB 


m: mean value, a 2 : variance, NSR: noise to signal ratio. SER: signal to error ratio. 


Table 4.4: Numerical results of 200 Monte Carlo runs for Example 4 
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Figure 4.11: Bode diagram of hi(s) and h 2 (s) in Example 4 
see that h 2 (s) has not only a noticeably lower frequency bandwidth, but also a smaller 
steady state gain. The results of 200 Monte Carlo runs are listed in Table 4.4, and 
the normalized bias (NB) and standard deviation (NSTD) computed from the table 
are plotted in Figures 4.12 and 4.13 respectively. In this case, the correspondence of 
NSRi and NSR 2 is 


relationship between NSRi and NSR 2 

noise 

level 1 

level 2 

level 3 

level 4 

NSRi 

0% 

5.03% 

10.10% 

20.04% 

nsr 2 

0% 

17.30%1 

34.72% 

68.85% 


Due to the huge amount of equivalent noise imposed on h 2 (s), both NB and NSTD 
values for h 2 (s) are approximately doubled comparing to Example 1. Like the results 
in Example 3 for system M 5 ), the NB numbers are almost diminished to be biasless, 
and the NSTD values are reduced to only half of those in Example 1. These results 
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Figure 4.12: Normalized bias vs. additive noise level NSRi 
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Figure 4.13: Normalized standard deviation vs. additive noise level NSRi 
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indicate the sharpest contrast between h\(s) and h 2 (s) among the four examples 


considered. 

Example 5: As the final example, we consider the system 


y( s ) 


hi{s)ui(s) + h 2 (s)u 2 (s) 


s 2 -f 3s + 8 


«i(«) + 


1.4 


s 2 + 2s + 


yU 2 (s). 


(4.28) 


The Bode diagrams of k^s) and h 2 {s) are plotted in Figure 4.11 from which we 



Figure 4.14: Bode diagram of hi(s) and h 2 (s) in Example 5 


can see that h 2 (s) has not only a noticeably lower frequency bandwidth, but also a 
higher steady state gain. The results of 200 Monte Carlo runs are listed in Table 4.5, 
and the normalized bias (NB) and standard deviation (NSTD) computed from the 
table are plotted in Figures 4.15 and 4.16 respectively. In this case, the corresponding 
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Figure 4.15: Normalized bias vs. additive noise level NSRi 
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Figure 4.16: Normalized standard deviation vs. additive noise level NSRi 
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hi(s) J 

1,1 ' 

h 2 (s) 



nisER 

true 0 

5 

3 

8 

1.4 

2 

1 


C*SER 

m 

4.9943 

3.0016 

7.9971 

1.3995 

1.99% 

0.99% 

0 % 

58.99dB 

CP 

0.00001 

0.00001 

0.00001 

0.00000 

0.00001 

0.00001 

0 % 

2.297dB 

m 

4.9908 

3.0013 

7.9893 

1.3937 

1.9887 

0.9954 

5 . 11 % 

44.76dB 

CP 

0.00226 

0.00097 

0.00141 

0.00045 

0.00156 

0.00022 

4.63% 

3.725dB 

m 

4.9816 

3.0012 

7.9708 

1.3760 

1.9572 

0.9821 

10.53% 

38.24dB 

CP 

0.00824 1 

0.00355 

0.00580 

0.00183 

0.00637 

0.00099 

9.54% 

3.682dB 

m 

4.9389 

2.9919 

7.8626 

1.2983 

1.8174 

0.9247 


30.85dB 

CP 

0.04197 

0.01485 

0.04599 

0.00973 

0.03006 

0.00521 

J 17.33% 

5.112dB 


m: mean value, o J : variance, NSR: noise to signal ratio, SER: signal to error ratio. 


Table 4.5: Numerical results of 200 Monte Carlo runs for Example 5 


relationship between NSRi and NSR2 is 


relationship between NSR X and NSR 2 

noise 

level 1 

level 2 

level 3 

level 4 

NSRi 

0% 

5.11% 

10.53% 

19.12% 

nsr 2 

0% 

4.63% 

9.54% 

17.33% 


From above form, we can see that additive noise levels in both channels are in the 
same range, which is analogous to Example 1. However, narrowing the bandwidth 
of h 2 (s) has made the NB and NSTD numbers for /t x (s) much smaller comparing to 
Example 1. For /i 2 (s), the NSTD numbers are also decreased and the NB numbers 
remain basically unchanged with respect to Example 1. Again, the iteration steps 
required for convergence of the AWLS/MFT algorithm in this case are noticeably 
fewer than those cases in Example 1 and 2. 
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4.3.4 Cross Comparisons and Comments 


In our arrangement, system hi(s) and its I/O pairs (ui(<), yi(t)) have been kept the 
same in all five examples; therefore, noise generated by equation (4.21) guarantees 
that the three additive noise levels (about5%, 10%, 20%) are physically in the same 
range in all five examples. After the five examples regarding different SISO subsystem 
combinations have been presented, some cross comparisons among them will be made 
here with the hope that further revelation could be made about the influence of 
different subsystem combinations on estimate bias and deviations. 

When there is no additive noise present in the output y(t), the parameters have been 
perfectly retrieved in all five MISO systems, no matter what combination it has. 

A narrower bandwidth as in Examples 3 and 5 has surprisingly improved both the 
speed and the quality of estimation on both h\(s) and h 2 (s) (especially for system 
/ii(s)). This has not been an anticipated result, because our modulating frequency 
range has been chosen based on the higher frequency bandwidth, so that the sys- 
tem h 2 (s), which has a narrower bandwidth, has been overmodulated with respect 
to its bandwidth. Pearson and Lee [39] had demonstrated that for the LS/MFT al- 
gorithm, this kind of overmodulating could incur estimation error. Therefore, only 
the contrary was expected at first. These “puzzling” results have been simulated 
many times and it has been found that the results are very much repeatable. This 
fact, on the other hand, is helpful in making the following two conclusions: (1) The 
AWLS/MFT algorithm is not as sensitive to the chosen modulating bandwidth or to 
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overmodulating as is the LS/MFT scheme"*; (2) The difference in bandwidth might 
actually make it easier for the AWLS/MFT algorithm to separate or distinguish the 
contributions from the two systems hi(s) and h 2 (s) to the single available output, 
and this might be more or less like a two-category classification problem in a pattern 
recognition framework in which the greater the distance between two clusters in the 
feature vector space, the better or more accurate the classification will be. 

When a lower gain h 2 (s) is used as in Example 2, both the NB and NSTD values 
for h 2 (s) were up slightly with respect to Example 1, but meanwhile the NB and 
NSTD numbers for hi[s) came down a little bit, correspondingly. From Example 3 
to Example 4, the gain of h 2 (s) was reduced from 0.6 to 0.3; here the adverse impact 
on the higher gain subsystem h\(s ) has not been noticed at all relative to that in 
Example 3, while the NB and NSTD numbers for h 2 {s) have been doubled. Therefore, 
lowering the gain of a SISO subsystem in the MISO identification problem will lower 
the accuracy of estimation for that subsystem, mainly because of the increase of its 
NSR 2 numbers. This is also consistent with Example 5, where increasing the gain of 
system h 2 (s) has improved its NSTD numbers slightly. 

Comparing to Example 1, the composite effect of narrowing the bandwidth and low- 
ering the gain of h 2 (s) relative to h- l (s) has been demonstrated more revealingly in 
Examples 4 and 5. In Example 4, the NB and NSTD values for h^s) have been 

decreased to half those in Example 1, and meanwhile these two numbers have been 
3 This property of AWLS/MFT will be further illustrated in the estimation of the Longitudinal 
dynamics of an F-18 aircraft (see Section 5.4.2). 
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almost doubled for h 2 (s). In Example 5, narrowing the bandwidth and increasing the 
gain of system h 2 (s) not only reduced the NB and NSTD numbers for hi(s), but also 
lower the NSTD numbers for system h 2 (s). 

Another measure of quality of parameter identification is the time domain perfor- 
mance of the estimated systems or models. The SER number mentioned in Sec- 
tion 4.3.2 is one such quantity characterizing the time domain performance of the 
model. In our Monte Carlo simulation studies, the mean and standard deviation 
of SER have been recorded and presented in Tables 4.1~ 4.5 from which the SER 
numbers in all five examples are fairly close to each other at each noise level NSRi; 
consistent with the NB and NSTD measures, Example 3 has about a 2dB edge over the 
equal bandwidth combinations. In order to provide some visual perception about the 
meaning of SER number, one typical Monte Carlo run in Example 1 with NSRi = 20% 
has been plotted in Figure 4.17. 


4.4 Concluding Remarks 

Without constraints among the denominator polynomials of a MIMO system, the 
parameter identification can be decomposed into several independent MISO sub- 
problems due to the joint cost function (4.15) stemming from the joint likelihood 
function of regression errors. The AWLS/MFT has been the core in forming the 
overall algorithm or procedure to solving this MISO identification problem. If the 
original form (4.1) is unmodulatable, it has to be converted into the modulatable 
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Figure 4.17: A typical realization of y(t), y(t), and e(t) = y(t) - y from Example 1 

form first and then apply the AWLS/MFT algorithm to estimate this higher order 
modulatable system. After that, as a model reduction tool, the AWLS/MFT scheme 
with specified lower orders can then be utilized for each subsystem to acquire the 
parameter identification of the original unmodulatable forms. Due to the fact that 
model reduction is carried out based on the estimate of the higher order system, it 
follows that the accuracy of the first AWLS/MFT estimate would be very crucial 
to the rest of the estimation. With the setup in our simulation studies, the results 
obtained in this higher order model stage are approximately a maximum likelihood 
estimate. 

The most intriguing phenomenon is the influence of the bandwidths of the two sub- 
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systems on the final quality of estimation. It is somehow contrary to the conventional 
wisdom to conclude that this difference of bandwidths actually benefits the estimate 
for both subsystems, even though the noise level NSR2 due to the narrowing band- 
width of subsystem /i 2 (s) has been almost doubled. Some underlying explanation 
from MFT itself or other approaches like artificial neural networks (ANN) should be 
targeted in the future research. Another implication of this phenomenon is that the 
AWLS /MFT algorithm is not sensitive to the modulating frequency bandwidth. 
Although the subsystem with lower gain has been seen as a drawback in our examples, 
the good time domain performance SER values have been persuasive enough to ensure 
its usability. Overall, the AWLS/MFT algorithm has been successfully applied to 
MISO parameter identification problems. 
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Chapter 5 


System Identification of the 
Longitudinal and Lateral 
Dynamics of an F-18 Aircraft 
Using the AWLS/MFT Algorithm 


5.1 Introduction 

5.1.1 Flight Variables Used in This Project 

Referring to Figure 5.1, the body axis system of an aircraft is taken with the center 
at the center of gravity (C.G.) of the airplane, OX forward, OY out the right wing, 
and OZ downward as seen by the pilot [4]. Most aircraft are symmetric with respect 
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X t 


* 

Figure 5.1: Sketch of earth axes and aircraft body axes 

to the OXZ plane. In order to describe the motion of the aircraft with respect to the 

earth or inertial space, a set of Euler angles, specifying the orientation of OXY Z with 

respect to an earth axis system OXeYeZe with its origin at the center of gravity of 

the aircraft and nonrotating with respect to the earth, can be used for this purpose 

(see Figure 5.1). We denote the three airplane body axis angular velocities as follows: 

Q : body axis pitch rate (rad/sec) 

P : body axis roll rate (rad/sec) 

R : body axis yaw rate (rad/sec). 

The transformation between the Euler angles ('P,0,<I>) and the angular velocities of 
the aircraft body axis (Q,P, R) can be written as [4] 

P = $ - 4> sin 0 
Q = © cos $ + 4 r cos 0 sin $ 
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R = — 0sin$ + \Hcos0cos$ 


Another two output quantities defined in Fig. 5.2 are (1) a : the angle of attack (rad), 
and (2) ft : the sideslip angle (rad) of the aircraft. In addition, three major control 

variables (driving forces) in the F-18 are 

Sh : symmetric horizontal tail deflection (rad) 

6 a : asymmetric aileron deflection (rad) 

S T : symmetric rudder deflection (rad) 

As mentioned earlier, we just need to be concerned about small perturbations, denoted 

as (a(t),ft(t),q{t),p(t),r(t),<t>{t)), around the equilibrium (or trim) operating point 

of flight, {a 0 ,fto,Qo,Po,R*,$oy- These small disturbances will be bounded with a 

set of linearized differential equations under a series of assumptions (see page 25 in 
1 As suggested by Dr. V. Klein, all the trim conditions in our studies will be determined by 

averaging the time domain sequences. 
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[4]); those coefficients (stability derivatives) of the differential equations will be the 
targets of this chapter. 


5.1.2 Criterion of Model Performance 


Far different from the simulation studies where the true parameters can be used luxu- 
riously to evaluate the quality of estimations by forming bias and standard deviations, 
when the original system is not known beforehand, the only seemingly reasonable al- 
ternative left is to compare the time domain performance of estimated models with 
the physical output. Traditionally, as used in the simulation studies, the signal to 
error ratio (SER), defined as 


SER = 20 • log j 0 


MQh ) 
Ili/W-s/WlhJ 


(5.1) 


could serve the purpose of evaluating the time domain performance, where y{t ) is the 
physical output instead of the “ideal” output (as in the simulation cases) and y(t) 
is the estimated output. Due to the long transient time of the aircraft system, the 
nonzero initial condition used in the MATLAB simulation routine LSIM() might play 
an important role in determining y(t) after the parameters are estimated, although 
this unknown initial condition has no impact on the quality of the estimations at all. 
In order to be as objective as possible to estimate y(t), a Luenberger Observer based 
initial condition (I.C.) estimation scheme was suggested by Pearson (see Appendix C) 
and it will be used as a standard tool in this chapter. 

As requested by NASA, the algorithms developed in this chapter will be tested by 


128 



estimating the parameters of theoretical models provided by NASA through Monte 
Carlo simulations at the suggested additive noise level. In these simulation cases, the 
performance criteria used in Chapter 4 for the estimated models will be employed 
without further ado. 


5.2 Identification of Longitudinal Dynamics 

5.2.1 Longitudinal Dynamical System 

As mentioned earlier, the longitudinal dynamics involve the response movement of 
the aircraft in the symmetric plane OXZ and its control signal. In this case, the 
input control signal would be the symmetric horizontal tail deflection 6h(t) an d the 
responses would be (i) the angle of attack Q'(t), and (ii) the body axis pitch rate <j{t) 
as shown in the block diagram (see Figure 5.3). Hence this is a single-input-and-two- 



(A(p), B,(p)) 

angle of attack a(t) 



horizontal tail deflection 



8h(t> 

(A(p),B e (p)) 

body axis pitch rate q(t) 




Figure 5.3: Block diagram of the longitudinal dynamics 
output (denote as SI20) system. The general dynamical differential equations are of 
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the forms 2 : 


A\ (p)a(t) = Bi(p)6h(t) (5.2) 

MpM*) = B 2 (p)6 h (t) (5.3) 

Mp) = Mp) = Mp) (5-4) 

where the orders of A(p), B\(p) and B 2 (p) are expected to be 2, 1 and 1 respectively; 
more importantly, there exists the physical pole constraint A\(p) = A 2 (p) = A(p) 3 , 
which make the algorithms developed in Chapter 4 inapplicable here. 

5.2.2 AWLS Algorithm for A Constrained SIMO System 

From the flight data received from NASA, the additive noise to 6h(t) is basically invis- 
ible, while the a(t ) and q(t) are contaminated with certain degrees of noise. Applying 
2 Basically this derives from 

x(t) = Ax(t ) + Bu(t), y(t) - Cx(t) 

where 

y(0 = ( ). u(t) = s h (t). 

According to 

y(s) = C{sI-A)~ 1 Bu(s), 

then 

det (si - A)y(s) = C(sl - j4) _1 det(s/ - j4)«(s). 

3 Private communication with Dr. V. Klein. 
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the measurement noise signal model and assuming that the additive noises to a(t) 
and q(t) are mutually independent white sequences, the joint likelihood cost func- 
tion, similar to (4.11), can be constructed. Therefore, for the constrained differential 
dynamical equation set: 

f A(p)a(t) = Bi(p)6h{t) (5 5 } 

\ A{p)q{t) = B 2 {p)6 h {t) 

consider the joint cost function: 

j(e u e 2 ) = (u - r 1 e 1 ) T w l -'(Y i - r 1 e 1 ) + u{y 2 - t 2 9 2 ) t w^(y 2 - v 2 e 2 ) (5.6) 


where constant v is a scaling parameter, 


0 i = 



and 9 a , 9 bl and 9 b2 are the parameter vectors comprising the coefficients of A(p), Bi(p) 
and B 2 (p) respectively. Partition the I\ and ^ according to: 


r\ = [r 01 ,r tl ], r 2 = [r a2 ,r b2 ] ; 


(5.7) 


such that equation (5.6) can be re-written as 


J(9aAM = (Yi + T a A-TbAyWr 1 (Y 1 + T ai 9 a -nA>)+ 
v{Y 2 + V a2 e a - T b2 9 b2 ) T W^(Y 2 + T a A - TbAi) 


The necessary conditions for minimizing J(9 a ,6 bl ,9 b2 ): 



5 


dj _ n ju_ _ 

ee hl “ u? dQ ~ 


lead to the following highly coupled equation set for 


+ vT T a3 W 2 -'T a2 )9 a = -r T a wr 1 Y 1 -vr T a2 w 2 -'Y 2 + 
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r«i*'f 1 r»A + (5.9) 


rjH7*r,A = rju’r'r, + rJwf'r.A (s.io) 

= rjwr'n + rJiVf'r.A (5,ii) 

Under the assumptions: 

• no measurement noise in the input signal 8h(t) 

• the additive noises on a(t ) and q(t ) are independently Gaussian distributed 
white processes with the variance ratios embedded in the scaling parameter v 

together with the pole constraint A/p) = A 2 (p), the weighting matrices Wi and W 2 
are identical, i.e., W\ = W 2 = W, so that updating and inverting W need be computed 
just once in each iteration. The scaling parameter v can be used to accommodate 
the measurement noise difference between the two output channels and to suppress 
or enhance the importance of the second output in the cost function. The following 
AWLS/MFT based algorithm has been employed to solve (5.9) rv (5.11): 

Algorithm 6 (Constrained SIMO AWLS/MFT Algorithm) 4 

1. Pick a scaling parameter value for v, v > 0, and estimate the initial value for 
0 a through the SISO system model: 

A(p)[o(t) + q(t)} = [Bi (p) + B 2 (p)]8 h {t) 

using the AWLS/MFT or WLS/MFT algorithm (see Section 2.2.2 and 2.2.3). 

4 Please refer to Section 2.3 for the detailed arrangement of real and imaginary quantities, espe- 
cially W, Wr and Wj in Equations (2.104) and (2.105), when the AWLS/MFT is implemented. 
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2. Compute the weighting matrix W and its inverse (see Section 2.2.2, 2.2.3 

and 2.3). 

3. Substitute the values for the pair {6 a ,W} into (5.10) and (5.11) and solve for 
the pair {9^,9^}. 

4- Estimate a new 0 a from (5.9) using the values for { 9 b 1 ,9)> 2 ,W} from the previous 

step. 

5. Check if the parameter value for the new 9 a has changed or not , based on a 
percent change in norm. If yes, go back to step 2, otherwise stop. 

6. Check the system output-signal-to-output error ratios S/E for the two models 
in (5.5) to see if they are in rough agreement with each other. If not, try a new 
value for v and repeat steps 1 ~ 6. 

5.2.3 Results Using Physical Flight Data 

Setup of Running AWLS/MFT Scheme 

With the sampling rate F s = 50Hz, the longitudinal maneuver flight data totalled 18 
seconds (or 900 signal points for each input and output signal). The other running 
parameters for the SIMO AWLS/MFT Algorithm 6 could be listed as 

• resolving frequency: /o = F s /900 = 0.056Hz 

• modulating bandwidth: Ft, = 0.5Hz 
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• highest modulating index: M = integer(Fb/F 0 ) = 9. 


As mentioned before, all the DC values of the flight data are subtracted as the trim 
operating conditions. Therefore, the time averages of all the data blocks (including 
both longitudinal and lateral channels) displayed in this Chapter are zero. 


Identified Longitudinal Dynamical Models 


As suggested by a theoretical model of NASA (see appendix A), the order of the 
model should be n = 2. However, two order assumptions, n = 2 and n = 3, were 
implemented for comparison purposes and the numerical values of the two estimated 
models through 20 AWLS/MFT iterations axe presented in equations (5.12) /~s_/ (5.15): 


a(s) _ -1.4769s -0.1361 

6 h (s) ~ s 2 + 0.7129s + 0.1804 
q(s) __ 0.1265s - 1.2313 

S h (s) ~ s 2 + 0.7129s + 0.1804 


(5.12) 

(5.13) 


n = 3 : 


<»(*) 

Sh(s) 

g(g) 

h(s) 


-1.0381s 2 - 0.5799s + 0.1013 
s 3 + 0.6487s 2 + 0.5507s + 0.0021 
-0.0995s 2 - 1.1270s - 0.5700 
s 3 + 0.6487s 2 + 0.5507s + 0.0021 


(5.14) 

(5.15) 


Their time domain performances are plotted in Figures 5.4 and 5.5 respectively. For 
n = 2 and v = 0.1, the estimated model has SER a = 7.89dB, SER, = 9.51dB. For 
n = 3 and v = 1, estimated model has performed even better, especially for channel 
cx(t). In these two cases, 20 iteration steps have been used. The evolution of SER 
numbers relative to the iteration steps during the process of Algorithm 6 can be 
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Figure 5.4: Time domain performance of estimated models, n 




time(s) 


Figure 5.5: Time domain performance of estimated models, n 
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further illustrated through the following discussions and Figure 5.6 where the case of 
n — 2 is plotted. 

The Influence of Scaling Factor v 

The overall impact of choosing different */, originally introduced to accommodate the 
variance difference of the two additive output noises, can be revealed partly from the 
SER-vs-iteration curves as typified in Figure 5.6, where a series of v values and their 

corresponding curves are displayed. At two extremes, e.g., from v — 10“ 3 — ► 10” 5 

12 

10 

S' 
v 

w 8 
OS 
w 

C/D 

6 
4 

0 5 10 15 20 25 

Iteration Steps 

Figure 5.6: The influence of u on the longitudinal identification 

and v = 10 3 — » 10 5 , the changes of the SER curves axe almost invisible as if they have 
reached two performance bounds. In order to enhance the performance of the a(t) 
channel, a smaller v is needed (equivalent to weighting the q(t) channel less in the joint 
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cost function (5.6)). Not surprisingly, reducing v will degrade the performance of the 
q(t) channel. Another factor demonstrating the affect of v is that in using v = 0.001, 
only 5 iteration steps are required to reach SER g = 9.74dB and SER a = 7.69dB 
(though SER a and SER„ are still varying rapidly), while 20 iterations are necessary 
when v = 0.1 is used. The subtlety of the choice in v seems far more profound than we 
have anticipated. The SER trade-off between the o:(f) and q(t) channels has been the 
dictating factor in deciding the v value in this case. However, this lack of elegancy in 
determining v has been far from devastating on account of the general insensitivity to 
i/, a quick finding for the v which balances the importance of the two output channels 
is almost trivial. Exploiting this flexibility will further enable the user to refine the 
results based on the particular physical setting. 

5.2.4 Simulation Results 

A numerical simulation study of Algorithm 6 was required by NASA, especially at a 
set of suggested additive noise levels. Using the theoretical model {A. 2) of NASA as 
the true system and driving it with the flight input data (horizontal tail deflection 
6 h {t))i 50 Monte Carlo runs are simulated for each of several additive noise levels. 
Like the previous simulation studies, the white Gaussian noise generator RANDNQ 
is used as the tool to manipulate the noise source. The additive noise levels (standard 
deviations cr’s) suggested by NASA for longitudinal quantities are 

• a q = 0.1 (degree) 
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• cr, = 0.05(degree/sec) 


Although it is a little bit ironic, o a = O.l(rad) and <7, = 0.05(rad/sec), instead of 
a a — 0.1 (degree) and o q = 0.05(degree/sec), were mistaken initially as the standard 
deviations of the additive noises, which had virtually amplified the original values 
with a scale 2n. In this case, the additive NSR numbers on a(£) and q(t ) had reached 
an appalling level, 200% and 95% respectively. Fortunately, the performance of Al- 
gorithm 6 is still quite promising (see Table 5.1, Figure 5.8 and 5.9). Realizing this 
embarrassingly erroneous scaling, it was decided that the noise level, a Q = 0.1 (degree) 
and a q = 0.05(degree/sec), would be used as a normalizing scale to introduce a “nor- 
malized” noise-to-signal ratio (NNSR), so that performance in a broader breadth of 
additive noise levels (at which the relative additive noise between a{t) and q(t) re- 
mains unchanged) could be seen (Figure 5.7). Under this normalizing configuration, 
the values suggested by NASA would be NNSR = 100% (its regular NSR a = 32% and 
NSR, = 15%). One typical realization of the Monte Carlo runs at NNSR = 100% is 
presented in Figure 5.8; it is seen that Algorithm 6 can provide nearly perfect results 
at the NASA suggested noise levels. Even at the devastating extreme NNSR = 628% 
where NSR a = 200% and NSR, = 95%, very impressive SER numbers still can be 
observed from Figures 5.7 and 5.9. 
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50 Monte Carlo Runs for Longitudinal Dynamics using constrained AWLS/MFT algorithm. Weight 0.25 


numerator 1 I numerator 2 | denominator | a | q 


true par I -0.0520 -1.4860 9 -1.5358 | -0.2435 


-0.2407 


-0.0568 -1.4842 


0.01168 0.01279 


-1.5333 


* I -0.0458 


0.01901 


-1.4862 -1.5355 


0.01843 I 0.01127 


0.00657 0.01006 I 0.00716 | 0.00362 


-0.2478 H 0.5111 


0.3579 | S/E(dB) I S/E(dB) 


0.3576 I 39.14 44.92 


3.091 3.482 


j|22Q£)| 


0.3592 


0.01094 0.00629 


-0.0531 

-1.4898 

0.02610 

0.03085 

-0.0698 

-1.4777 

0.05507 

0.06405 

-0.0454 

-1.5089 

0.09989 

0.13374 


-1.5356 I -0.2497 


-1.5313 I -0.2328 


0.5128 1 0.3595 


0.02062 0.00904 


0.03754 I 0.05004 I 0.03557 I 0.01733 


-1.5214 


0.06864 I 0.03403 


NNSR is normalized NSR by the case o„ = 0.1 (deg) |NSR=32%) and o, = 0.05(deg/sec) {NSR=15%}. 
•means the conditions suggested by NASA. 



Table 5.1: Numerical results of 50 Monte Carlo runs for the longitudinal dynamics 


pitch rate 



angle of attack 
NASA suggested case 


200 300 400 500 600 700 

NNSR (%) 


Figure 5.7: SER numbers at different additive noise levels 
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Figure 5.8: Time domain performance of Algorithm 6 at NNSR = 100% 
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Pitch Rate (rad/s) Angle of Attack (rad) 


aa = 0.l(rad) ideal output 

NNSR = 628%,NSR=200% noisy output 

0-4 -S/E=1 9.08dB . ; • estimated output 



Time(s) 


Figure 5.9: Time domain performance of Algorithm 6 at NNSR = 628% 
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5.3 Identification of Lateral Dynamics 


5.3.1 Lateral Dynamical Systems 




(A(p), B/p), BJp)) 

sideslip angle P(t) 

aileron deflection 5 a (t) 






(A(P), B zi (p). BJp)) 


body axis roll rate p(t) 







(A(p), B 3i (p), BJp)) 


body axis yaw rate r(t) 

rudder deflection fy(t) 






(A(p), B 4i (p), BJp)) 


Euler roll angle 4> (t) 



r 


Figure 5.10: Block diagram of the lateral dynamics 

The longitudinal dynamics of aircraft seem less complex than the lateral dynam- 
ics, which involve more inputs and outputs. As suggested by NASA, three 5 major 
responses are (1) the sideslip angle /?(f); (2) the body axis roll rate p(t)’, and (3) 
the body axis yaw rate r(t ) as plotted in the block diagram (see Figure 5.10). The 

5 In an early communication with NASA, they named these three quantities (/?(<), p(t), r(t)) as 
the outputs and our research was carried out under this configuration. Later on, from the theoretical 
model sent by NASA, the Euler roll angle <j>(t) was found to be embedded in the model as well, but 
it is loosely coupled with the other quantities. In this section, all the analysis will be presented in a 
two-input-and-four-output (2140) framework, from which the 2130 scheme can be derived trivially. 
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mathematical differential models describing them can be expressed as 


{ Ai(p)m = B n (p)Sa(t) + B u (p)6 r (t) 

I M (p)p(<) = B 2 l(p)S a (t) + B 2 2 {p) 8 r(t) (5 16 \ 

1 A 3 (p)r(<) = S 3 ,(p)*.(i) + £32(p)*r(i) 

[ Aj(p)</>(t) = 7?4l(p)<5o(0 + B42{p)6 T {t) 

The orders of {B i3 {p)]i = 1,2, 3,4; j = 1,2} and {A,(p);i = 1,2, 3, 4} are expected 
to be 3 and 4 respectively. The attachment of the physical pole constraint A x (p) = 
A 2 (p) = A 3 (p) = A 4 (p) — A(p) to the above differential equation set poses a genuine 


constrained MIMO system identification problem. 


5.3.2 AWLS Algorithm for A Constrained MIMO System 


Similar to the longitudinal case, the additive noises to both control signals (8 a (t),8 r {t)) 
are visually nil. Resorting to the measurement noise signal model and assuming the 
additive noises to (/?(<), p(t), r(<), <f>(t)) are mutually independent white noises, the 
joint likelihood cost function like (4.10) of Chapter 4 can be easily formed. Hence, 


for the constrained dynamical system: 

( A(p)m = B n (p)6 a (t) + B 12 (p)8 r (t) 

A(p)p{t) = B 2 \{p)8 a {t) + B 22 (p)8 r (t) 

A(p)r(t) = B 3 i(p)8 a (t) + B 32 (p)8 r (t) 

, A{p)(j>(t) = B 4 i(p)8 a (t) + B 42 (p)8 T (t) 

consider the joint cost function: 

+vi(y 2 - y 2 o 2 ) t w 2 \y 2 - r 2 0 2 ) 

+u 2 (y 3 - r 3 e 3 ) T w^(Y 3 - r 3 0 3 ) 

+v 3 (y 4 - y 4 o 4 ) t w 4 \y 4 - r 4 0 4 ) 

where v x , v 2 and v 3 are the scaling factors. The 0* vectors can be defined by 


0i = 
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where 9 a , 9 bl , 9 b2 , 9b 3 and 9b, are the parameter vectors comprising the coefficients 
of A(p ) and {Bn(p),Bi 2 ]i = 1,2, 3, 4} respectively. Partition the IT, IT, T 3 and T 4 
according to 


IT = [r 0l ,r 6l ] ; r 2 = [r aa ,r 6a ] ; r 3 = [r a3 ,r fc3 ] ; r 4 = [r fl4 ,ru 


such that equation (5.18) can be re-written as 

./(0.AAAA) = (n + r ai 0 o -r 6l ^ 1 ) r ip 1 - 1 (r 1 + r ai 0 a -r 6l 0 6l ) 

+»i(Y 2 + r a2 9 a - T h 9 b2 ) T W 2 - 1 (Y 2 + T a2 9 a - T h 6 h ) 
+v 2 (y 3 + r a3 0 a - ita)^ 1 ^ + r as 9 a - r b3 9 b3 ) 
+u 3 (Y 4 + r a ,9 a - r b ,9 b ,) T W 4 -'(Y< + T a ,9 a - Y b ,9 b ,) 

(5.19) 


The necessary conditions for minimizing J(9 a ,9 bl ,9 b2 ,9b 3 ,9b,): 


M. _ o 

d$a ’ 


d J n d J a dJ a dJ a 

d6 hl - dd b2 ~ dO b 3 “ u ’ d$ bi - u 


lead to the following highly coupled equation set for {9 a ,9 bl ,9 b2 ,9 b3 ,9b,}: 


(r£wf 1 r « 1 + Vl T T ai w^T a2 + V2 y lw^r a3 + = 

-T^Wr'Yi - I 'xTl 2 W^Y 2 - u 2 T T a3 W^Y 3 - v^W^Y, 


+rr, Wr'YbAi + viTlW 2 -'Y b2 9b 2 + v 2 Y T W^ 


'n n i _ _ rT 




r£wr'r 4 A = rlW{'Y,+rlwr'r.,e. (5.21) 

TlWj'r h e h = rlw.;'Y 2 + riw.;'r„0, (5.22) 

rjH-T'riA = rj^-'n + rj^-'r.,#. (5.23) 

rJw.r-nA = rju'.-'n + rj^-T.,^ (5.24) 


Under the assumptions: 
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• no measurement noises in the input signals 6 a (t ) and 6 r (t) 


• the additive noises on /3(t), p(t ), q(t) and <f>(t) are independent Gaussian dis- 
tributed white processes with the variance ratios embedded in the scaling pa- 
rameters V\ , v 2 and v 3 

and taking account of the pole constraint A/p) = A 2 (p) = /^(p) = A 4 (p), the 
weighting matrices Wi, W 2 , W 3 and W 4 are identical, i.e., Wi = W 2 — W 3 = W 4 — W, 
so that updating and inverting W need be computed just once in each iteration. 
The scaling parameters 1/1 , v 2 and 1/3 can be used to accommodate the intensity in 
measurement noise differences among the four output channels and to suppress or 
enhance the importance of the p(t), r(t ) and <f>[i) outputs in the cost function. The 
following AWLS/MFT based algorithm has been employed to solve (5.20) ~ (5.23): 

Algorithm 7 (Constrained MIMO AWLS/MFT Algorithm) 6 

1. Pick nonnegative scaling parameter values for (r'l, v 2 , V3), and estimate the ini- 
tial 0 a through the MISO system model: 

A(p)[0(t) + p(t) + r(t) + <f>{t)\ = [Bn(p) + B 2 /p) A B 3 /p) + £? 4 i(p)] 6 a (t)-|- 

[■Si 2 (p) + B 22 (p) + B 32 (p) + 5 42 (p)]^ r (0 

using AWLS/MFT or WLS/MFT algorithm (see Section 2.2.2 and 2.2.3). 

2. Compute the weighting matrix W and its inverse (see Section 2.2.2, 2.2.3 

and 2.3). 

6 Please refer to Section 2.3 for the detailed arrangement of real and imaginary quantities, espe- 
cially W , Wr and Wi in Equations (2.104) and (2.105), when the AWLS/MFT is implemented. 
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3. Substitute the values for the pair { 6 a , W } into (5.21) ~ (5.23) and solve for the 
parameter set {8b , , 0b 2 , 0b 3 , 9b t } • 

4. Estimate a new 8 a from equation (5.20) using the values for (^,^,^,^1 W] 
from the previous step. 

5. Check if the parameter value for the new 8 a has changed or not, based on a 
percent change in norm. If yes, go back to step 2, otherwise stop. 

6. Check the system output-signal-to-output-error ratio S/E for the three models 
in (5.17) to see if they are in rough agreement with each other. If not, try a 
new value for the triple (i>i, V 2 , 1/3) and repeat 1 ~ 6. 

5.3.3 Results Using Physical Flight Data 

Setup of Running the AWLS/MFT Scheme 

With the sampling rate F s = 50Hz, the lateral flight test data totalled 20.48 seconds 
long (or 1024 signal points for each I/O channel). The running parameters used for 
MIMO AWLS/MFT Algorithm 7 are listed below: 

• resolving frequency: f 0 = F s j 1024 = 0.0488Hz 

• modulating bandwidth: Fb = 0.5Hz 

• highest modulating index: M = integer(Fb/F 0 ) = 10. 
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Identified Lateral Dynamical Models 


In this case, three output channels (/?(<), p{t), r(f)) and two inputs (6 0 (i), 6 r (t)) are 


used to form a two-input-and-three-output (2130) system with constrained 4th order 


denominator polynomial. After 22 iterations as shown in Figure 5.11, Algorithm 7 


with V\ = 1 and V 2 — l 7 was truncated and resulted in the following estimated models: 


P( s ) 


p(s) 


r{s) 


0.2217s 3 - 1.1464s 2 + 0.5981s - 0.8248 f , , 
s 4 + 0.4633s 3 + 2.9945s 2 + 1.0374s + 1.5875 al5j 
-0.0544s 3 + 0.3066s 2 - 0.1994s + 0.3975 
s 4 + 0.4633s 3 + 2.9945s 2 + 1.0374s + 1.5875 * ri5j 
-4.2991s 3 + 0.8664s 2 - 7.8554s + 3.8868 
s 4 + 0.4633s 3 + 2.9945s 2 + 1.0374s + 1.5875 ’ 
0.4304s 3 - 0.9799s 2 + 0.7119s - 2.5624 ^ 

s 4 + 0.4633s 3 + 2.9945s 2 + 1.0374s + 1.5875 
-0.3512s 3 + 1.0683s 2 - 1.2603s + 2.5927 
s 4 + 0.4633s 3 + 2.9945s 2 + 1.0374s + 1.5875 ' 
-0.1687s 3 - 0.4299s 2 - 0.2223s - 1.1016 
s 4 + 0.4633s 3 + 2.9945s 2 + 1.0374s + 1.5875 ' 


(5.25) 


(5.26) 


(5.27) 


Approximately 3 minutes computation time was needed to calculate these mod- 
els using an IBM 486/33 machine. The poles of the above models are located at 
(-0.1312 ± 1.3899* ) and (-0.2235 ± 0.8010*). The time domain performance of the 
above models are plotted in Figure 5.12. As pointed out in the NASA model (Ap- 
pendix B), the fourth state (</>(t)) could be used as a fourth output. The identifica- 
tion for this configuration (equivalent to a 2140 system like (5.17) was implemented 
as well by using the Algorithm 7; notice that the resulting SER numbers for the 

three channels p(t), r(t )) slipped from (9.41, 11.62, 11.60)dB to the current 

7 Other v values have been tested and, like the longitudinal case, the results are not sensitive to 


the change of v - 
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Figure 5.11: SER numbers vs. iteration steps for the lateral dynamics 

(7.80, 6.50, 10. 44)dB with V\ = 1, Vi = 1 and Us = 1. However, it should not be too 
surprising to see that the estimated 2130 linear models axe capable of interpreting the 
physical data better than those identified with a theoretically linearized 2140 model 
structure. 


5.3.4 Simulation Results 

Similar to the longitudinal case, numerical simulation studies on Algorithm 7 were 
requested by NASA at a set of suggested additive noise levels. Using the theoretical 
model (B. 2) of NASA as the true system and driving it with the recorded flight input 
data (<5) a (t) and 6 r (t)), 5(3 Monte Carlo runs axe simulated for each additive noise level 
on this 2140 system. Like the previous simulation studies, the white Gaussian noise 
generator RANDN0 is used as the tool to manipulate the noise source. The additive 
noise levels (standard deviation <r’s) suggested by NASA for the lateral dynamics are 
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• crp = 0.1 (degree) 


• <t p = 0.05(degree/sec) 

• a T = 0.05(degree/sec) 

• <7^, = 0.01 (degree). 

Due to the initial misusage of unit (rad), we actually have tested Algorithm 7 at 
other larger additive noise settings. In order to represent those results as well, the 
NASA suggested noise level, the crp — O.l(degree), a p = 0.05(degree/sec), a r = 
0.05(degree/sec) and a $ = 0.01 (degree), will be used as a normalizing scale to form 
the normalized noise-to-signal ratio (NNSR). Then, the NASA suggested noise level 
would be equivalent to NNSR = 100% (NSR^ = 51%, NSRp = 8%, NSRr = 14%, 
and NSRtf, = 1%). As mentioned earlier, the Euler angle was excluded from 
our initial studies on a constrained 2130 system. Here in the simulation study, the 
underlying system is known to be the linear 2140 system (B. 2). The AWLS/MFT 
algorithms for both the 2130 and 2140 models have been tested for this ideal 2140 
system; their resulting numerical values are listed in Tables 5.2~ 5.3. As plotted in 
Figure 5.13, except for the the yaw rate r(t) at the low noise side, the time domain 
performances of the two models are fairly close to each other. One typical realization 
for the 2140 model at the NASA suggested noise level is plotted in Figure 5.14, 
which confirms that the identified model from Algorithm 7 does give a very good 
time domain performance. 
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0.4978 


0.4957 


0.03881 


denominator 


2.0426 0.5494 


2.0449 0.5458 


0.02845 0.04582 


-0.0049 


-0.0051 


0.00373 


Signal-to-Error Ratio (dB) 


P 


20.58 28.60 28.41 26.72 


4.479 4.438 4.331 4.233 


NSR* 

NSRp 

NSRr 

NSR+ 


5% 63% 

0 . 6 % 


0.4875 I 2.0450 I 0.5398 I -0.0048 | 15.81 24.29 | 24.48 25.05 


0.05198 I 0.04271 I 0.07106 I 0.00490 I 4.840 4.052 4.319 | 6.126 


0.4884 

2.0434 

0.5414 

-0.0049 

0.09070 

0.06864 


0.00874 

0.4658 

2.0590 

0.5233 

-0.0040 

0.18730 

0.13720 

0.25562 

0.01793 

0.4475 

2.5177 

0.6517 

-0.0126 

0.24609 1 

0.47250 

0.51007 

0.04366 


12.79 

20.58 

20.94 

19.88 

4.503 

4.234 

4.556 

5.900 

6.44 

12.80 

13.03 

12.74 

6.137 

5.228 

4.861 

6.458 

-1.28 

1.78 

7.31 

5.76 

6.822 

7.378 

6.203 

7.422 


160% 

25% 

44% 

3.1% 

320% 

50% 

88% 

6.2% 


2140 


true 6b, 


mean 


std 


B„(p) 


0.0005 -1.4171 -0.3838 -0.0010 


-0.0019 -1.4074 -0.3905 0.0078 


0.05820 0.05482 0.14572 0.04695 


Bdp) 


0.3607 0.1013 


-0.0096 NNSR 


-0.0132 63% 


0.00963 I 0.01307 I 0.03070 I 0.02189 


-0.0045 


0.09684 


-1.4179 


0.08393 


-0.4042 


0.23247 


-0.0078 


0.08610 


0.01657 I 0.02239 0.05456 


-0.0107 | 100% 
0.03436 I 


0.0062 

-1.4258 

-0.3950 

-0.0160 

0.16831 

0.16294 

0.35702 

0.09945 

-0.0734 

-1.4869 

-0.5868 

-0.0257 

0.31107 

0.22358 

0.78659 

0.21128 

0.0652 

-1.4847 

0.0506 

0.1420 

0.67416 

0.64572 

1.63521 

0.54428 


0.0101 

0.3624 

0.1057 

-0.0085 

0.02840 

0.04245 

0.07059 

0.05433 

0.0242 

0.3630 

0.1394 

-0.0311 

0.05095 

0.06657 

0.17025 

0.10618 

0.0099 

0.2721 

0.0562 

-0.0262 

0.11964 

0.21488 

0.32295 

0.23633 


2140 


true 0^ 


mean 


std 


B*(p) 


-3.3714 -0.1695 -0.4805 0.0151 | 0.3206 


-3.3631 -0.1555 -0.4740 0.0153 | 0.3197 


0.06355 0.12431 0.15993 0.08470 I 0.01706 


B^Cp) 


-0.0933 -0.9442 


-0.0936 -0.9429 


0.02368 0.05818 


0.0297 I NNSR 


0.0330 


0.03427 


-3.3561 


0.09697 


-0.1425 


6269 


-0.4498 


0.25232 


- 0.0022 


0.11859 


0.3210 


0.02634 


-0.0956 


0.03352 


-0.9391 


0.08786 


0.0404 | 100% 
0.04286 I 


mean 

-3.3627 

-0.1348 

-0.4482 

0.0176 

0.3197 

-0.1004 

-0.9454 

0.0348 

157% 

std 

0.14111 

0.29050 

0.40911 

0.18397 

0.04014 

0.05711 

0.14626 

0.07324 


mean 

-3.3072 

-0.0475 

-0.4288 

-0.0950 

0.3184 

-0.1109 

-0.9042 

0.0703 

314% 

std 

0.30623 

0.5539 

0.90239 

0.38923 

0.09137 

0.10302 

0.31985 

0.14981 


mean 

-2.4245 

-0.1986 

0.5914 

-1.7 182^ 

0.1630 

0.0565 

-0.8759 

0.6018 

628% 

std 

1.22100 

0.96646 ^ 

3.09565 

1.68237 

0.22311 

0.31769 

0.65002 

0.61901 



Table 5.2: 50 Monte Carlo runs for the 2140 system, continued to next page 
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2140 B„(p) B ?2 (p) 

true 6,3 0.0512 | 0.0443 -0.2122 | -0.0528 -0.2480 | -0.1205 | -0.4240 | -0.1037 NNSR 

mean -0.0604 0.0467 -0.1875 -0.0496 -0.2491 -0.1190 j -0.4290 -0.1004 63% 

"ltd 0.03487 0.02961 0.07364 0,02448 0.00571 "001420" 0.01677 0.01306 

mean* | -0.0530 | 0.0604 | -0.2092 | -0.0472 | -0.2485 | -0.1217 | -0.4239 | -0.1027 | 100% 
std* 0.05347 0.03910 0.12307 0.0427lT 0.00804 ~0.01814 0.03026~ 0.02251 

mean | 0.0510 | 0.0624 | -0.1960 | -0.0442 | -0.2477 | -0.1220 | -0.4263 | -0.1019 | 157% 

"ltd 0.06528 0.08418 0.17229 0.05617 0.01150 0.03605 0.04080 0.03247 

mean 0.0668 0.0445 j -0.2110 -0.0420 -0.2512 -0.1116 -0.4306 -0.1010 314% 

“ltd 0.15520 0.13645 0.43128 0.12453 0.02849 0.05995 0.10341 0.06591 

mean 0.0202 0.04505 -0.2021 0.0045 j -0.2403 -0.1176 -0.5743 -0.1330 628% 

"ltd 0.29986 0.25990 0.78443 0.31560 0.05071 0.09594 0.19534 0.16327 


2140 B 4 ,(p) B 42 (p) 

true 6 m 0.0000 | -3.3567 | -0.1568 [ -0.5413 0.0000 | 0.2495 | -0.1278 | -1.0658 NNSR 

mean 0.0168 -3.3495 j -0.0751 -0.5449 -0.0050 0.2492 -0.1485 -1.0566 63% 

“ std 0.11747 0.05168 0.65505 0.05851 ~OQ4259~ 0.01555 0.18734 0.09533 

mean* I 0.0290 I -3.3485 I 0.0128 I -0.5366 | -0.0100 I 0.2498 I -0.1747 I -1.0471 | 100% 

std* 0.18092 0.08935 0.97541 0.07951 0.06590 0.02154 028038 0.14951 

mean | 0.0223 | -3.3494 | -0.0326 | -0.5339 | -0.0078 | 0.2490 | -0.1630 | -1.0532 | 157% 

std 0.29690 0.13044 1.62360 0.14707 0.10710 0.03651 0.46487 0.24128 

mean j 0.0866 j -3.3152 j -0.3207 | -0.5324 j -0.0288 j 0.2403 j -02543 j -1.0216 j 314% 

“std 0.66787 0.33048 3.57656 0.29508 0.24097 0.07876 1.02168 0.54732 

mean 0.4973 -3.3082 j 1.9556 -0.8647 j -0.0669 0.0767 -0.3828 -1.1494 [ 628% 

“std 1.46643 0.52769 6.73534 1.52819 0.43575~ 0.35374 ' 1,73141 1.30480 

NNSR is normalized NSR by the case: 

o„ = 0.1(deg){NSR=320%}, a t = O.Q5(deg/sec) {NSR=50%}, 
a R = 0.05(deg/sec) {NSR=90%} and a* = 0.01(deg) {NSR=4%} 

‘means die conditions suggested by NASA. 

Table 5.2: 50 Monte Carlo runs for the 2140 system. 
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denominator 


0.4978 2.0426 0.5494 | -0.0Q19 


0.5896 2.0898 0.7455 0.0365 


0.07444 0.03430 0.16395 0.02451 


Signal-to-Error Ratio (dB) 


P P R I 9 


20.24 24.16 16.50 


5.568 4.139 


NSRp NSRp 
NSRr I NNSR 


32% I 5% I 63% 


0.6384 I 2.1117 I 0.8310 I 0.0609 


2285 0.08427 0.28007 0M930 


17.10 23.01 14.28 


4.372 6.131 3.322 


51% 8% 100% 
14% 


0.6435 


0.12532 


0.8666 


0.08223 I 0.26781 


1.1456 


2130 


true 0 bi 


mean 


std 


0.7178 

2.3007 

0.3840 

0.22805 

0.6107 

3.7172 

0.27017 

0.97554 


j j> 

| 0.0005 

-1.4171 


0.0645 


0.04421 


0.11557 


0.3034 


0.23535 


14.06 I 19.42 I 13.52 


4.570 4.561 3.136 


12.25 

10.42 

5.364 

3.237 


0.93 I 4.01 I 6.83 


3.953 4.083 3.021 


13% 157% 


160% 25% 314% 

44% 


320% 50% 628% 

88 % 


-o.ooto 


0.06130 


-0.3838 -0.0010 


-1.4150 I -0.5423 -0.0124 


0.06025 0.22370 0.06063 


0.0098 I 0.3594 I 0.1389 I -0.0086 


0.01000 0.01475 0.04175 0.02495 


-0.0073 


0.11714 


-1.4320 


0.10742 


-0.5980 


0.30615 


-0.0302 


0.08245 


0.3651 I 0.1529 I -0.0028 | 100% 


0.02092 I 0.03301 


0.0109 


0.17574 


-1.4299 


0.17090 


-0.5937 


0.42460 


-0.0530 


0.11126 


0.0119 I 0.3621 I 0.1620 0.0042 


0.03038 0.04210 0.07919 0.06003 


mean 

-0.0562 

-1.5255 

-0.8527 

0.0050 

0.0248 

0.3439 

0.2060 

0.0051 

314% 

std 

0.31627 

0.23401 

1.09784 

0.25488 

0.05642 

0.08442 



mean 

0.0432 

-0.9580 

0.2054 

-0.0341 

0.0459 

0.0591 

0.3657 

0.0812 

628% 

std 

0.61253 

0.52412 

1.99219 

0.82101 

0.13793 

0.16181 

0.58878 

0.27237 



2130 


true 0 b2 


mean 


std 


-3.3714 


-3.3101 


0.13358 


B^p) 


-0.1695 -0.4805 


-0.4506 I -0.2708 I -0.0814 


0.25478 0.58695 0.13645 


B a (p) 


0.3206 -0.0933 -0.9442 


0.3045 -0.0635 -1.0048 


0.03132 0.02744 0.14083 


0.0297 


-0.0203 


0.06415 


mean* I -3.3026 I -0.6240 I -0.3516 | -0.0944 \ 0.3061 


std* I 0.17993 0.38391 0.65234 0.25284 | 0.04237 


-0.0440 -1.0032 


0.04667 


-0.0695 


mean 

-3.2706 

-0.5983 

-0.1552 

-0.1597 

std 

0.24672 

0.41293 

0.95006 

0.30031 

mean 

-2.8878 

-0.7190 

0.5547 

-0.8582 

std 

0.68298 

1.29488 

2.20614 

0.98283 

mean 

-0.4345 

-0.6585 

3.0656 

-4.6326 

std 

1.88277 

1.08985 

5.19007 

2.24115 


0.2936 

-0.0503 

-1.0345 

-0.0493 

0.05082 

0.05359 

0.20757 

0.11916 

0.2331 

-0.0143 

-1.0607 

0.0725 

0.14829 

0.17885 

0.50853 

0.35707 

-0.0942 

0.4349 

-0.9247 

1.1387 

0.33222 

0.39915 

1.10963 

0.84094 


Table 5.3: 50 Monte Carlo runs for the 2130 system, continued to next page 
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| 2130 | 

B»(P) 

Bsi(p) 


1 Ena 

0.0512 

0.0443 

-0.2122 

-0.0528 

-0.2480 

-0.1205 

-0.4240 

-0.1037 

NNSR 

mean 

0.02974 

0.0407 

-0.3028 

-0.0124 

-0.2403 

-0.1358 

-0.4105 

-0.1658 

63% 

std 

0.04208 

0.03095 

0.10482 

0.04578 

0.00818 

0.01889 

0.02032 

0.04579 

mean* 

0.0059 

0.0414 

-0.3671 

0.0523 

-0.2341 

-0.1470 

-0.3993 

0.0523 

100% 

std* 

0.06223 

0.04826 

0.20007 

0.08745 

0.01266 

0.03102 

0.03453 

0.09147 

mean 

-0.0033 

0.05834 

-0.3993 

0.0616 

-0.2321 

-0.1517 

-0.3957 

-0.2117 

157% 

std 

0.07577 

0.08486 

0.20906 

0.0992 

0.01485 

0.03745 

0.03863 

0.08203 

mean 

-0.0302 

0.0937 

-0.5789 

0.1915 

-0.2275 

-0.1728 

-0.4083 

-0.3194 

314% 

std 

0.14836 

0.15328 

0.45579 

0.26600 

0.02724 

0.10844 

0.09409 

0.24733 

mean 

-0.4184 

0.3270 

-1.9588 

1.0482 

-0.1594 

-0.2433 

-0.5940 

1.0482 

628% 

std 

0.41884 

0.27371 

1.42094 

0.91631 

0.07261 

0.10470 

0.35619 

0.44713 


NNSR is normalized NSR by the case: 

Op = 0.1(deg){NSR=51%}, o P = 0.05(deg/sec) {NSR=8%}, 
o R = 0.05(deg/sec) {NSR=14%} and o 4 = 0.01 (deg) {NSR=1%J 


♦means the conditions suggested by NASA. 


Table 5.3: 50 Monte Carlo runs for the 2130 system. 


5.3.5 A Brief Comment on Minimal Realization 


As seen in the above simulation studies, the true system is a fourth order state space 
model, but Algorithm 7 only returns the I/O differential models instead of state-space- 
looking I/O model like (B.l). One question that automatically occurs to our mind 
is whether the estimated transfer function model, at the NASA suggested noise level 
through Algorithm 7, could be minimally realized using a fourth order state space 
model. This has been in no sense a trivial question to answer. Assuming distinct 
poles, the well-known Gilbert’s Diagonal Realization Scheme [15] first expands the 
transfer function matrix H(s) into partial fractions: 


H(s) 


NW 

d(s) 


e; n,s- n 

s — \i j s — A, 


(5.28) 
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where the denominator polynomial is 


d(s) = f](s - A,) 

1 

and {R} are the residue matrices. Further denote 

Pi = the rank of Ri 

and then Gilbert’s method says that the minimal realization has order 

n = J2pi (5.29) 

1 

So the question of obtaining a minimal realization of the transfer function model 
H(s) has boiled down to the determination of the residue matrices {it,} and their 
ranks. As suggested by Professor A.E. Pearson, employing the partial fraction ex- 
pansion algorithm in [19] and the rank determination scheme based on the singular 
value decomposition technique [10], some verifications on {Ri} and {p,} have been 
computed. Define the tolerance for rank determination 8 by: 

&SVD — ' |[it,||oo 

where d max is the precision of estimating the matrix elements, i.e., the maximum 
standard deviation value among all the entries comprising each TV, matrix. This is 
obtained from the std values in Table 5.2~ 5.2. At the NASA suggested noise level, 
our calculations on the models estimated by Algorithm 7 with different additive noise 
realizations show that all the {it,} have persistently been diagnosed as rank one 
8 All the singular values less than 6svd will be dropped. 
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matrices. This means that the estimated transfer function model can be realized with 
a fourth order state space model and demonstrates consistency with the theoretical 
state space model. 


5.4 Miscellany 

Inasmuch as many other issues pertaining to the F-18 dynamics have been investigated 
as well, some relatively-peripheral-but-intriguing results should be worth mentioning 
briefly in this section. Hopefully, some merit of the SISO AWLS/MFT Algorithm 2 
could be further divulged through handling the physical flight data. 

5.4.1 Identification of the Actuator System 

The F-18 horizontal tail deflection 8h(t) is activated by the longitudinal pilot stick 
movement rjh(t) through a mechanical actuator system. Using the SISO AWLS/MFT 
Algorithm 2 in Section 2.2.3, the dynamics of this actuator system can be modeled 
as a linear second order system resulting in the transfer function: 

6 h (s) _ —0.0494s 2 - 0.0462s - 0.0340 

rjhls) ~ s 2 + 0.2586s + 1.2987 ' ( ' 

The time domain performance of the above model has been plotted in Figure 5.15 
from which we can see that the model yields an impressive time domain fit. As to 
why a second order model is used, several other order model structures were also 
tested and tabulated. The model (5.30) results from using the parsimony principle: 
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Figure 5.15: Time domain performance of the estimated longitudinal actuator 

pick the lowest order model among all the models which give a reasonable time do- 
main performance. Together with the time domain performance, other approaches 
like statistically checking the whiteness of the residual process W -1 ^ 2 (Y — T6) using 
Lemma 1 with a 95% confidence bound [5] [52] had also been adopted to determine 
the model structures in our early research report to NASA. 


5.4.2 Sensitivity to the Chosen Modulating Bandwidth Fb 

As mentioned earlier, our experience is that the AWLS/MFT Algorithm is far less sen- 
sitive to the chosen modulating bandwidth Fb than either the LS/MFT or WLS/MFT 
algorithm. Pearson and Lee in [39] found that overmodulating (broader modulating 
bandwidth than the system bandwidth) could exaggerate the estimation errors when 
LS/MFT is used. In our initial studies on the longitudinal dynamics, a fourth or- 
der SISO model was once built to link the longitudinal pilot stick movement to the 
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body axis pitch rate q(t). The physical spectrum of the body pitch rate q(t) has a 
bandwidth around 0.5Hz. During this model-building process, a set of modulating 
bandwidths ranging from 0.3Hz~1.0Hz had been chosen for LS/MFT, WLS/MFT us- 
ing both the unconstrained and constrained AWLS/MFT algorithms. The results are 
summarized in Table 5.4. It is seen that the models for the LS/MFT and WLS/MFT 


Sensitivity of dil 

fferent algorithms to the modulating bandwidth Fb 

algorithm 

F b (Hz)/ 

/m 

AWLS 

constrained 

AWLS 

unconstrained 

WLS 

LS 

0.3/12 

■table 

ste=5.84dB 

■table 

s/e-5.84dB 

unstable 
s/e= -4.78dB 

■table 

s/e= -6.79dB 

0.4/16 

•table 

s/e -6. 83d B 

■table 

&/e=6.83dB 

unstable 
s/e- 2.09dB 

unstable 
s/e« -3.39dB 

0.5/20 

•table 

s/e=7.13dB 

■table 

s/e-7.13dB 

■table 

s/e>1.86dB 

unstable 

s/e-1.53dB 

0.6/24 

•table 

s/e=6.91dB 

unstable 

s/e-7.56dB 

■table 

S/e=3.08dB 

unstable 
s/e* “7.41 dB 

0.7/28 

table 

s/e-6.51dB 

unstable 

s/e-8.07dB 

•table 

s/e-3.32dB 

unstable 
s/e- -53.22dB 

0.8/32 

■table 

ste«6.41dB 

unstable 

s/e-8.l0dB 

unstable 
s/e« -2 7. 63d B 

unstable 
s/e« -41 04dB 

0.9/36 

■table 

s/e=6.56dB 

unstable 

s/e-7.94dB 

unstable 
s/e- -33.98dB 

unstable 
s/e- -27.60dB 

1.0/40 

■table 

s/e=6.47dB 

unstable 

s/e=8.00dB 

unstable 
s/e= -61 ,97dB 

unstable 
s/e= -121.1dB 


Using physical flight data to build a fourth order model linking the 
longitudinal pilot stick movement (input) and the body pitch rate (output). 


Table 5.4: Sensitivity of MFT algorithms to the chosen modulating bandwidth 


are either stable with poor SER numbers or unstable with unacceptable SER’s. On 
the other hand, the estimated systems from both the constrained and unconstrained 
AWLS /MFT algorithms yield much more consistent time domain performance with 
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no drastic fluctuations like the LS/MFT and WLS/MFT, no matter whether the esti- 
mated models are stable or not. All these help claim that the AWLS/MFT algorithm 
not only performs better, but also has less sensitivity to the pre-chosen modulating 
bandwidth Fb- 
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Appendix A 


Longitudinal Model from NASA 


For simulation purpose, the following linearized theoretical longitudinal models cou- 
pling a(t) and q(t ) are given by NASA: 


f a{t) = -0.1692a(t) +0.95609(/) -0.0520^(t) 

{ q(t) = -0.3141a(t) -0.3403q(t) - 1.5358S h (t) * 


The corresponding constrained differential dynamical model is 


where 


/ A(p)a(0 = 
l A(p)q(t) = B 2 (p)6 h (t) 


A{p ) = ( 1.0000 0.5095 0.3579 ) 



( B i(p)\ = ( -0.0520 -1.4860 \ ( P \ 

{ B 2 {p) ) \ -1.5358 -0.2435 ) ^ 1 ) ' 


(A. 2) 


(A.3) 


(A.4) 
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Appendix B 


Lateral Model from NASA 


For simulation purpose, the following linearized theoretical lateral dynamical models 


coupling /?(f), p{t), r(t) and <f>(t) are given by NASA: 


m ) 


-0.0432 0.4066 -0.9108 0.1003 \ 


m \ 

i>(t) 


-4.4588 -0.4441 0.2692 0 


Pit) 

r(t) 


0.2239 -0.0057 -0.0106 0 


r(t) 



0 1.000 0.2867 0 ) 


4>{t) / 


0.0005 0.0098 \ 

-3.3714 0.3206 / S a (t) \ 

0.0512 -0.2480 \ $ r (t) / ' 

0 0 / 


The corresponding constrained differential dynamical model is 

r A(p)P(t) = B u (p)8 a (t) + B u (p)6 T (t) 
A(p)p(t) = B 2 i{p)8 a {t) + B 22 (p)8 T (t) 

A(p)r(t) = B 31 (p)6 a (t) + B 3 2(p)8 r {t) 
k A(p)<f>(t) = B 4 i(p)6 a (t) + B 4 2(p)8r(t) 


where 


A(p) = ( 1.0000 0.4978 2.0426 0.5494 


-0.0049 ) 


p 4 \ 

P 3 

P 2 

P 

\1 / 


(B.l) 


(B.2) 


(B.3) 
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Bn(p) > 


0.0005 -1.4171 -0.3838 -0.0010 > 


P 3 \ 

^2l(p) 


-3.3714 -0.1695 -0.4805 0.0151 


P 2 

B 3 i(p) 


0.0512 0.0443 -0.2122 -0.0528 


P 

V B 4 i(p) ) 


0.0000 -3.3567 -0.1568 -0.5413 > 


1 y 

b 12 (p) \ 


f 0.0098 0.3607 0.1013 -0.0096 \ 

/P 3 \ 

B 22 (p) 


0.3206 -0.0933 -0.9442 0.0297 


P 2 

B 32 (p) 


-0.2480 -0.1205 -0.4240 -0.1037 


P 

\ B 42 (p) / 


V 0.0000 0.2495 -0.1278 -1.0658 / 


V 1 / 


(B.4) 


(B.5) 
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Appendix C 


Initial Condition (I.C.) Estimate 
via Luenberger Observor 


The following algorithm for estimating I.C. is suggested by Professor A.E. Pearson. 
Given (A(p), B(p)) for the differential model A(p)y(t) = B(p)u(t), 0 < t < T, with the 
I/O data pair [u(t),t/(f)j on [0,T], we want to find x 0 = x(0) for the state realization 


x(t) = Ax(t) + Bu{i) 

y(t) = Cx(t) and xq = x(0) 


where (A, B, C) is any observable state space realization for 


|ifl 

A(*)’ 


i.e., 


C(sI-A)~ 1 B = 


B(») 

Ms) 


The I.C. estimate can be constructed through the following steps: 
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1. Subtract off the zero-state response y u (t) = Ce A ^ t ^Bu(t)(1t from y(t): 

y(t ) — y u {t) = Ce Ai x o I.C. response 

2. Define the reversed I.C. response 

*(*) = V{T ~ t) - y u (T - t) 

-- Ce~ At e AT x o 

= Cx R (t) (C.l) 

where x R (t ) = e~ At e AT x 0 satisfies: x R (t) = — xk(0) = e AT x 0 , and more 
importantly x R (T) = x 0 . Hence, we estimate x R (t) given z(t) via an observer 
as follows: 

x R (t) = - Ax R (t ) + L(z(t ) - Cx R {t )) 

= — {A + LC)x R (t) + Lz(t) 

with xh( 0) an arbitrary estimate. In order to carry out the error analysis, let 
x(t) = x R (t) — x R (t), then 

i(t) = -Ax R {t)-[-Ax R (t) + L{Cx R (t)-Cx R (t))) 

= —A(x R (t) - x R (t)) - LC(x R (t) - x R (t)) 

= ~{A + LC)x{t) 

which implies 

x(t) = e~ {A+LC)t x(0) — ► 0 
provided — {A + LC ) is Hurwitzian. 

3. Design the gain matrix L such that e ~( A+LC '> T « 0 and use x R (T) as an estimate 
of x 0 , i.e., x r (T) rs x 0 . 
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Appendix D 


Some Results in M.S. 
Differentiation and Integration 


The following results are cited mainly from the Chapter 4 of [53]. 

D.l Mean Square Continuity 

Definition 1 (continuous in mean square) A second-order stationary (or wide 
sense stationary) s.p. X(t), t € T , is continuous in mean square, or m.s. continuous, 
at a fixed t if 

l.i.m.r^ 0 X(t + T) = X(t) (D.l) 

Theorem 1 (continuity in mean square criterion) A second-order s.p. X{t), 
t € T, is m.s. continuous at t if, and only if, T(t,s) = £{A(f)A(s)} is continuous 
at (t,t). 
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Definition 2 (m.s. continuous on an interval) If a second-order s.p. X(t), t € 
T, is m.s. continuous at every t £ [^1,^2] C T, then X(t) is m.s. continuous on the 
interval [^i, <2] • 


D.2 Mean Square Differentiation 


Definition 3 (m.s. derivative) A second-order s.p. X(t), t E T, has a mean 
square derivative (or m.s. derivative) X(t) at t if 


X(t + r) — X(t) ^ 
Li.m.j^ o = X(t) 


(D.2) 


Definition 4 (differentiable on an interval) If a second-order s.p. X(t), t E T, 
is m.s. differentiable at every t E [<i, ^2] C T, then X(t) is m.s. differentiable on the 
interval [<1, <2] * 


Theorem 2 (criterion of m.s. differentiable) A wide sense stationary second- 
order process X(t), t E T, is m.s. differentiable if, and only if, the first- and second- 
order derivatives of the correlation function T(f,s) = ^{^(^^(s)} exist and are 
finite at t — s = 0. 


D.3 Properties of Mean Square Derivatives 

1. Mean square differentiability of X(t) at t E T implies m.s. continuity of X(t) 
at t. 
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2. The m.s. derivative X(t) of X(t) at t € T, if it exists, is unique. 

3. If ,X(£) and Y(t) are m.s. differentiable at t € T, then the m.s. derivative of 
a.X’(t) + bY(t) exists at t and 

4 W) + fer(t)] = aX(t) + bY(t) (D.3) 

at 

where a and b are constants. 

4. If an ordinary function f(t) is differentiable at t <E T and X(t ) is m.s. differen- 
tiable at t G T, then f(t)X(t) is m.s. differentiable at t and 

>W01 = ®*W + /M^ (D.4) 

D.4 Mean Square Riemann Integration 

Consider a collection of all finite partitions {p n } of an interval [a,b]. The partition 
p n is defined by the subdivision points tk, k = 0, 1,2, such that 

a — 1 ^ 1 2 ^ ... ^ — b 


Let 


A n = max(fi - 

k 

and let t' k be an arbitrary point in the interval i , ^ fc)* Let X(t) be a second-order 

s.p. defined on [a, b] C T. Let f(t,u) be an ordinary function defined on the same 
interval for t £ [a, 6] and Riemann integrable for every u £ U. We form the random 
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variable 


Y.( u) = £/«,“)*«)(<»-**-•) 

k = 1 

Since L 2 — space is linear, Y^(u) is an element of the L 2 — space. It is a random variable 
defined for each partition p n and for each u £ U. 

Definition 5 (m.s. Riemann integral) If, for u eU, 

U.m.n^oo^^oY^u) = Y(u) 

exists for some sequence of subdivisions p n , the s.p. Y(u), u € U , is called the mean 
square Riemann integral, or m.s. Riemann integral of f(t,u)X(t), over the interval 
[a, 6], and it is denoted by 

Y(u) — I* f(t,u)X(t)dt (D.5) 

Ja 

It is independent of the sequence of subdivisions as well as the positions of the t' k € 

[f *-!,<*)• 

Theorem 3 (integration in mean square criterion) The s.p. Y(u), u € U, de- 
fined by Eq. {D.5) exists if and only if the ordinary double Riemann integral 

f f f{t,u)f(s,u)Txx(t,s)dtds 
Ja Ja 

exists and is finite. 

D.4.1 Properties of mean square Riemann integrals 

1. Mean square continuity of X(t) on [a, b] implies m.s. Riemann integrability of 
X{t ) on [a, 6]. 
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2. The m.s. integral of A"(t) on an interval [a, 6], if it exists, is unique. 


3. If X(t) is m.s. continuous on [a, 6], then 

I f b X(t)dt\< f \\X{t)\\dt <M(b-a) 

\\Ja Ja 


where 

te[a } b\ 

4. If the m.s. integrals of X(t) and Y(t) exist on [a,c], then 

[ C [aX{t) + 0Y(t)]dt = a f X(t)dt + 0 f Y(t)dt 

Ja Ja Ja 

[ C X(t)dt = f X(t)dt + r X(t)dt, a < b < c 

Ja Ja Ja 

5. If X(t) is m.s. continuous on [a,t] C T, then 

V(t)= f* X(t)dt 

Ja 

is m.s. continuous on T; it is also m.s. differentiable on T with 


(D.6) 


(D.7) 

(D.8) 


m = m 


Corollary 1 (Leibniz Rule) If X(t) is m.s. integrable on T and if the original 
function f[t , s) is continuous on T xT with a finite first partial derivative df(t, s)jdt, 
then the m.s. derivative of 

Y{t)= f f(t,s)X(s)ds 
J a 

exists at all t £ T, and 

m= ( d . 9 ) 
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Corollary 2 (Integration by Parts) Let X(t) be m.s. differentiable on T and 
let the ordinary function /(/,s) be continuous on T x T whose partial derivative 
df(t,s)/dt exists. If 

Y(t)= f f(t,s)X(s)ds (D.10) 

J a 

then 

Y(t) = f(t , s)X{s) |‘ - j (‘ ?IfA X ( s )ds (D.ll) 

Let f(t,s) = 1 in Eqs. (D.10) and (D.ll); we have 

X(t)-X(a)= f X(s)ds, M€T 

J a 

This property is seen to be the m.s. counterpart of the fundamental theorem of the 
ordinary calculus. 

D.4.2 Means and Correlation Functions of M.S. Riemann 
Integrals 

Corollary 3 (means and correlations of m.s. Riemann Intetral) If the m.s. 
Riemann integral 

K(u)= [ b f(t,u)X(t)dt 

Ja 

exist, then 

E{Y(u)} = f b f(t , «)£{*(*)}<** (D.12) 

• Ja 

and 

D{y(u)y(u)}= f f f(t,u)f(s,v)E{X(t)X(s)}dtds (D.13) 

Ja Ja 
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